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1.0  INTRODUCTION 


The  recent  profound  shift  in  the  global  balance  of  power  in  favor  of  the  United  States  of  America 
has  had  major  repercussions  on  Strate^c  Defense  Initiative  (SDI)  planning.  In  particular,  the 
focus  has  shifted  from  the  provision  of  protection  for  the  United  States  against  a  massive  raid, 
involving  possibly  thousands  of  reentiy  vehicles,  to  defense  against  a  much  more  limited  attack 
which  could  now,  however,  be  launched  from  any  part  of  the  world.  Additionally,  the  United 
States  is  seeking  to  protect  its  forces  and  allies  overseas,  and  in  the  task  of  missile  detection  and 
tracking,  allowance  must  now  be  made  for  trajectories  which  can  begin  and  end  in  almost  any 
inhabited  area  of  the  globe.  Thus  SDI  demands  on  surveillance  technology  have  been  significantly 
expanded  [1]. 

Space-based  imaging  systems  will  play  a  vital  role  in  the  surveillance  task.  In  the  detection, 
discrimination  and  subsequent  tracking  of  hostile  objects,  image  quality  in  particular  will  be  of 
crucial  importance.  Considerations  of  size,  weight,  and  cost  will,  however,  impose  strict  limits  on 
the  unaided  performance  of  the  optical  hardware.  In  addition,  the  sensor's  operational  capabilities 
may  be  impaired  by  ^stem  aberrations  and  degradations  which  may  be  inherent  (as  in  the  case  of 
the  Hubble  telescope),  induced  by  the  stresses  of  launch  or  thermal  fluctuations  in  orbit,  or  simply 
the  result  of  aging  processes  in  the  space  environment.  If  these  defects  can  be  properly 
characterized,  however,  newly  developed  algorithnu  can  be  used  to  compensate  for  them  and  thus 
restore  the  image.  The  resolving  power  of  the  system  can  also  be  extended  in  this  way. 

In  this  report  a  description  is  first  given  of  a  typical  scenario.  The  potential  imaging  problems  are 
then  examined,  the  mathematical  background  is  discussed,  and  the  innovative  algorithms  which 
have  been  developed  for  correcting  and  enhancing  the  performance  of  the  imaging  sensor  are 
described.  Some  results  are  included  from  current  simulations  based  on  the  parameters  of  the 
typical  scenario.  The  complexity  involved  in  the  calculations  is  assessed,  together  with  the 
associated  storage  requirements.  Attention  is  drawn  to  wider  applications  of  the  techniques 
developed  and  knowledge  gained  in  the  course  of  this  work,  and  finally  recommendations  are 
made  concerning  the  design  attributes  of  an  orbiting  sensor  network  appropriate  to  the  SDIO 
missiotL 
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2.0  THE  SCENARIO  AND  THE  PROBLEM 


2.1  The  Defensive  Task 

Orbiting  rings  of  small,  Ughtweight  surveillance  satellites  (the  so-called  'brilliant  eyes')  provide  the 
surveillance  function  from  high  earth  orbit.  They  communicate  among  themselves,  or  via  a 
command  and  control  system  to  produce  a  multi-sensor  map  of  the  target  area.  The  number, 
nature  and  trajectories  of  objects  which  are  deemed  to  be  potential  threats,  m  the  atmosphere  or  in 
space,  are  identified.  Compact  guided  missiles  (the  so-called  'brilliant  pebbles')  are  on  stand-by 
and  are  dispatched  by  the  command  and  control  system  to  intercept  specific  objects  among  the 
group  of  threats.  The  missiles  destroy  or  damage  their  targets  through  high-speed  impact. 

The  optical  requirements  for  the  two  groups  of  satellites  (the  'eyes'  and  the  'pebbles')  differ  as  a 
consequence  of  their  intended  functions.  A  pebble  missile  will  be  initially  directed  towards  an 
intended  target,  which  may  be  moving  at  a  very  high  relative  speed,  and  must  be  capable  of 
acquiring  and  locking-on  to  it  as  rapidly  as  possible.  Moreover,  the  missile  must  carry  out  these 
operations  in  a  potentially  very  cluttered  environment,  and  be  capable  of  making  critical  final 
corrections  to  ensure  success.  Thus  simple  quadrant  detectors  and  centroid  techniques  would  not 
be  appropriate  except  very  near  the  target.  The  compleaty  of  the  pebble's  sensor  system  can  be 
greatly  reduced  by  designing  it  to  receive  initial  orientation  instructions  from  the  control  center, 
then  to  achieve  lock-on  and  make  course  corrections  autonomously.  The  eyes,  on  the  other  hand, 
must  be  capable  of  surveying  a  wide  field  of  view  and  of  achieving  high-resolution  imaging  to 
provide  the  information  for  detection,  identification  and  tracking.  The  computational  techniques 
discussed  in  this  document  are  designed  to  enhance  images  of  the  type  obtained  by  an  eye,  and 
would  be  applicable  also  to  a  pebble's  optical  system,  if  on-board  resources  permitted. 


2.2  Typical  Diffraction-Limited  Image  Properties 

On  the  basis  of  certain  assumptions  about  the  characteristics  of  the  optical  system  and  of  the 
target,  we  can  make  predictions  of  its  performance.  We  shall  estimate  system  resolution 
capabilities,  fields  of  view  and  available  image  acquisition  time  for  a  moving  target,  and  photon 
flux  into  the  image. 
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We  make  the  following  assumptions  for  system  parameters: 


Wavelength  X 

Aperture  d 

Focal  length  f 

Pixel  width  p 

Range  to  target  z 


10  pm 
250  mm 

1250  mm  (F/5  system) 
25  pm 
2000  km 


The  range  quoted  is  a  maximum  for  an  oibhal  height  of 700  km. 
Then: 


Rayleigh  resolution  length  at  image  plane 

1.22  Xfi'd  =  61  pm 

Resolvable  separation  at  range  z 

= 

1.22  Xz/d  =  98  m 

Airy  disc  diameter 

s 

2.44  Xfid  =  122  pm 

Pbcels  across  Airy  disc  = 

2.44  Xfidp  »  5 

Pixel  angular  field  of  view  (FOV) 

p/f  =  20prad 

FOV  at  2000  km 

sa 

40  m 

For  a  1000-element  linear  staring  array: 

Total  angular  FOV  =  20  mrad  x  20  prad 

Linear  FOV  at  2000  km  =  40kmx40m 

In  order  to  calculate  the  maximum  time  available  for  acquisition  of  a  single  image,  we  need  an 
estimate  of  the  maximum  velocity  of  the  object  across  the  field  of  view.  For  re-entrant  ballistic 
trajectories,  we  shall  assume  that  this  is  of  the  ordor  of  the  orbital  speed  for  a  low  earth  orbit,  and 
use  a  value  of  7  km  s*^ .  Then  the  image  of  a  target  at  2000  km  range  crosses  one  pbcel,  having  a 
40  m  FOV,  in  approximately  6  ms.  and  the  resolvable  distance  of  98  m  in  14  ms.  A 
reasonable  image  acquisition  time  might  be  half  of  this;  i.e.,  7  ms. 
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We  now  calculate  the  photon  flux  at  the  detector  for  a  typical  target. 


Planck's  formula  for  the  power  radiated  per  unit  area,  Aly,  within  the  frequency  interval  Av,  by 
a  black  body  at  temperature  T  is  [2] 


AT  -  o^Av 

”  c*  {exp(ho  /  kX)}  - 1 

where  h  is  Planck's  constant  (6.6  x  10*^  Is),  k  is  Boltzman's  constant  (1.38  x  10-^  JK'>)  and 
c  is  the  velocity  of  light  (3  x  10*  ms-^). 

For  X,  =  10  |im,  u  =  c/X  =  3  x  10‘3  s-‘. 

Suppose  the  body  is  radiating  at  300  K  and  that  the  radiation  is  detected  over  a  bandwidth  of  9- 
llfim.  Then  the  corresponding  frequency  increment,  Au,  is  O.OxlO^^s*^  and  from  Planck's 
formula,  the  power  per  unit  area  AI^  is  approximately  63  W  m-^. 

Suppose  further  the  target  is  radiating  into  4n  steradians  from  an  area  of  1  m*  .  Then  the 
radiated  power  into  the  250  mm  aperture  is  approximately  6  x  10-^^  W .  The  energy  of  a  photon 
at  10  pm  is  2  X  10-**  J .  Hence  the  photon  flux  through  the  aperture  is  3  x  10*  s-^ .  For  an 
acquisition  time  of  7  ms ,  the  number  of  photons  contributing  to  the  image  will  be  21  x  10* .  If 
the  target  is  essentially  a  point  source,  these  photons  are  distributed  over  the  point  spread 
function  of  the  system.  About  86%  of  this  energy  goes  into  the  central  lobe,  the  Airy  disc,  which 
covers  about  25  pixels.  Hence  the  average  energy  per  pbcel  in  the  Airy  disc  is  about  720  photons. 
The  center  pixel  would  in  fact  contain  about  3800  photons. 

Poisson  statistics  govern  the  energy  distribution  in  the  image  around  its  expected  value,  and  the 
tignal-to-noise  ratio  (SNR)  is  given  by  the  square  root  of  the  mean  signal  value  in  photons.  For 
the  average  pixel  energy  of  720  photons,  the  corresponding  SNR  is  27;  for  the  center  pbcel  its 
value  is  about  61. 
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2.3  Deficiencies  in  the  Sensor  System 


In  the  absence  of  any  other  degrading  effects,  the  performance  of  an  optical  system  is  restricted 
ultimately  by  the  effects  of  diffraction.  The  finite  extent  of  the  exit  pupil  imposes  a  fundamental 
upper  limit  on  the  system's  spatial-fi-equency  response.  It  is  unlikely,  however,  that  the  image 
quality  of  an  operational  space-based  system  will  approach  the  theoretical  limit  very  closely  at  any 
time.  It  is  possible  that  the  design  or  construction  will  be  flawed,  as  in  the  case  of  the  Hubble 
telescope,  through  defective  manufacture,  assembly  or  quality  assurance  procedures.  In  addition, 
it  is  highly  likely  that  the  sensor  will  become  degraded  during  the  mission  lifetime;  faults  may  be 
induced  by  the  stresses  of  launch  and  deployment,  and  aging  of  the  equipment  (caused,  for 
example,  by  radiation  or  thermal  effects)  will  almost  certainly  compromise  component 
performance.  The  detector  itself  will  impose  limitations;  for  example,  where  a  CCD  array  is  used, 
information  is  lost  in  the  inter-pixel  areas,  and  image  energy  is  integrated  over  the  active  area  of 
each  pbcel.  Other  degrading  fiictors  will  include  defective  pbcels,  noise  in  the  CCD  array  and 
electronic  subsystems,  and  the  possibly  obtrusive  contributions  of  the  earth-space  background. 
Success  in  detection,  identification  and  subsequent  tracking  will  depend  critically  on  the  levels  of 
noise  and  clutter  in  the  images,  and  robustness  is  of  fundamental  importance  in  any  image- 
processing  scheme.  The  methods  of  image  restoration  considered  here  were  originally  aimed  at 
achieving  performance  beyond  the  conventional  diffraction  limit  [3,  4],  but  are  in  fact  capable  of 
compensating  simultaneously  or  separately  for  aberrations  induced  by  the  optical  components  and 
for  the  limitations  of  the  detector.  They  were  designed  to  be  robust  and  also  possess  valuable 
noise-suppression  properties. 
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3.0  THE  SOLUTION 


3.1  Mathematical  Bacl^ound 

The  initial  assumption  is  made  that  the  effect  of  the  optics  can  be  described  as  a  possibly  time-  and 
shift-variant  blurring  of  the  image  due  to  diffraction  and  aberrations.  Thus,  in  general,  the  point 
spread  function  vnll  change  across  the  sensor  field  of  view.  It  can  be  assumed,  however,  that  at 
any  given  time  the  point  spread  function  is  effectively  constant,  i.e.,  shift-invariant,  over  some 
localized  area,  and  then  undergoes  a  discrete  change  into  neighboring  areas.  This  is  the  case  with 
the  Hubble  telescope.  It  will  also  be  assumed  that  the  set  of  point  spread  flmctions  is  known  or 
can  be  measured.  Under  these  circumstances  the  point  spread  function  is  said  to  be  locally  shift- 
invariant,  and  the  image  is  created  by  the  summation  across  the  entire  field  of  view  of  the  set  of 
localized  point  spread  functions  convolved  with  objects  in  the  corresponding  parts  of  the  field. 
Mathematically,  a  convolution  can  be  written  as  a  Fredholm  integral  equation  of  the  first  kind. 
The  solution  of  this  class  of  equation  is  known  to  be  ill-posed  and  numerically  difficult  [S,  6]. 
Additionally,  it  is  anticipated  that  the  image  will  be  spatially  sampled  by  a  solid-state  sensor  which 
introduces  spatial  integration,  discretization  and  associated  noise  processes.  Thus,  the  Fredholm 
integral  sum  representing  the  continuous  image  can  be  rewritten  as  a  matrix  equation,  to  which 
modem  rignal  processing  techniques  derived  fi’om  advanced  linear  algebra  can  be  applied.  In 
general,  the  presence  of  the  sensor  noise  takes  the  measured  image  out  of  the  span  of  the  columns 
of  the  kernel  matrix,  which  is  typically  highly  ill-conditioned;  additional  techniques  derived  fi’om 
regularization  theory  are  required  to  restore  stability  to  the  reconstmction.  By  introducing  a 
suitable  error  criterion,  images  can  be  constmcted  which  are,  in  terms  of  the  chosen  criterion, 
closer  to  the  undistorted  geometrical  image  of  the  object  than  the  detected  image  data. 

The  point  spread  function  of  a  weU  designed  optical  system  is  normally  approximately  invariant 
over  large  segments  of  the  image  plane.  Note  that  ima^ng  performance  (and  therefore  the  point 
spread  fimcdon)  commonly  changes  over  relatively  wide  fields  of  view;  for  example,  the 
resolution  of  typical  photographic  lenses  is  maricedly  less  near  the  edge  of  the  field  than  near  the 
center.  However,  over  some  limited  region  the  point  spread  function  will  be  essentially  shift 
invariant,  although  this  is  not  mathematically  necessary,  and  one  may  speak  of  an  effective  point 
spread  function  over  this  re^on  with  a  transfer-function  band-limit  determined  by  the  exit  pupil  of 
the  qrstem.  System  aberrations  and  pbcel  integrarion  serve  to  reduce  the  system  response  within 
this  passband. 
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To  the  extent  permitted  by  the  noise  in  the  image,  these  in-band  eflfects  are  relatively 
straightforward  to  remove  for  shift  invariant  systems;  a  typical  approach  would  be  to  employ  a 
Wiener-type  inverse  filter.  However,  detector  pbcellation  and  the  finite  aperture  of  any  system  set 
resolution  limits  not  easily  overcome,  and  a  method  for  achieving  spectral  extrapolation  has  to  be 
devised. 

The  spatial  spectrum  of  the  object  is  the  Fourier  transform  of  its  amplitude,  in  the  coherent  case, 
or  its  intensity,  in  the  incoherent  case.  If  the  object  is  known  to  be  of  finite  extent,  its  Fourier 
transform  is  an  analytic  function,  and  the  out-of-band  part  of  the  spectrum  can  in  principle  be  fully 
recovered  by  analytic  continuation  [7]  of  the  image  spectrum  after  removal  of  any  in-band 
distortion.  The  inverse  Fourier  transform  of  this  extended  spectrum  would  then  yield  a  perfect 
image  of  the  ori^al  object.  Equivalently,  one  could  attempt  to  solve  directly  the  equation 
describing  the  ima^g  process.  This,  however,  involves  the  inversion  of  an  ill-conditioned 
matrix,  and  the  restoration  process  is  intrinsically  unstable,  even  small  amounts  of  noise  rendering 
the  results  meaningless.  These  difSculties  may  be  surmounted  by  applying  the  methods  of 
regularization  theory  [8],  developed  to  deal  with  ill-posed  problems  of  this  type,  and  the 
procedures  which  have  evolved  in  the  course  of  this  SDIO/ISTO  program  are  based  on 
construed  least-squares  methods  in  which  a  regularization  parameter  plays  an  essential  role. 
Stability  in  the  restored  image,  which  is  computed  via  the  regularized  pseudoinverse  of  the 
ima^ng  matrix  (see  Appendix  A),  is  controlled  by  this  parameter.  Its  optimal  value  depends  on 
the  signal-to-noise  ratio  in  the  data.  If  it  is  not  possible  to  select  this  optimum  value  firom  prior 
knowledge  of  the  system  characteristics  and  the  object  scene,  several  techniques  are  available  for 
its  estimation  firom  the  data  themselves  [9,  10].  Faster  computation  of  the  estimate  is  possible  if 
the  singular  value  decomposition  of  the  imaging  matrix  is  already  available. 

An  efBcient  algorithm  should  exploit  known  properties  of  the  target  as  much  as  possible.  In  the 
midcourse  phase  the  targets  of  interest  radiate  thermally  at  approximately  300  K  and  the  images 
are  real,  nonnegative  distributions.  The  re-entry  vehicles  are  relatively  small,  presenting  an  area 
o^  say,  0.2S  n^,  and  are  spinning,  whereas  the  bus  is  larger,  with  an  area  of  several  square  meters. 
All  targets  are  in  ballistic  trajectories.  The  extent  of  the  targets  is  much  smaller  at  most  ranges 
than  the  resolution  limit  of  the  optics  described  in  Section  2.2,  so  that  they  typically  act  as  point 
sources. 
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Effective  image  restoration  requires  as  good  an  estimate  as  possible  of  the  imaging  system 
characteristics,  careful  design  of  the  jrithm,  appropriate  stabilization  of  the  algorithm  against 
noise,  and  high  computational  efficiency.  The  algorithms  which  have  been  developed,  and  these 
which  are  still  under  development,  differ  primarily  in  the  extent  to  which  the  target  characteristics 
can  be  identified  and  incorporated  into  the  algorithm,  the  dominant  noise  mechanisms  in  the 
imaging  process,  and  the  degree  to  which  the  computational  efficiency  can  be  optimized. 

The  most  straightforward  approach  to  obtaining  the  necessary  imaging  system  characteristics  for 
image  enhancement  would  be  to  perform  pre-launch  measurements  on  a  test  rig.  However, 
allowance  should  be  made  for  the  possibility  of  launch  misalignments  and  in-orbit  changes 
occurring;  see  Section  2.3.  This  can  be  done  in  two  ways.  The  first  is  to  modify  the  algorithnu 
to  allow  for  uncertainty  in  the  calibration  data  as  well  as  noise  in  the  image  data;  the  algorithms 
then  depend  on  'total*  least-squares  methods  [11].  The  second  is  to  test  the  point  spread  function 
periodically  and  update  the  system  data.  In  this  context  it  may  be  noted  that,  if  the  eyes  are 
capable  of  active  surveillance,  it  may  be  possible  to  use  the  light  source  in  one  to  calibrate 
another. 

3.2  Image  Restoration  Algorithms 

Reconstruction  algorithms  can  be  divided  into  two  major  classes;  those  for  which  the  object  can 
be  allowed  to  take  both  positive  and  negative  values,  and  those  for  which  non-negativity  is  an 
important  constraint.  Algorithms  in  which  non-negativity  is  imposed  on  the  reconstruction  are 
inevitably  iterative  and  require  greater  computational  effort  to  perform  the  minimization.  The 
lower  computational  burden  of  the  non-iterative  form  of  the  algorithms  is  available  if  non¬ 
negativity  does  not  have  to  be  imposed;  this  may  be  the  case,  for  example,  for  relatively  bright 
objects  which  stand  out  clearly  fi-om  their  surroundings.  For  a  group  of  small  isolated  objects,  on 
the  other  hand,  use  of  the  non-negativity  criterion  can  be  a  powerful  means  of  improving  object 
location.  Thus  the  methods,  which  also  provide  smoothed  numerically  stable  solutions  to  the 
matrix  equation,  are  well-suited  for  SDI  space-based  surveillance  applications.  Algorithms  have 
been  constructed  for  a  number  of  specific  purposes,  including  the  correction  of  optical 
aberrations,  for  achieving  resolution  beyond  the  diffiaction  limit  of  the  optical  system  and  for  the 
recovery  of  apparently  lost  sub-pbcel  detail.  The  performance  of  some  of  these  algorithms  will  be 
illustrated  with  simulated  and  laboratory  images.  (See  also  Appendbc  B.)  Preliminary  studies 
have  also  been  made  of  image  data  acquired  during  space  launches  by  the  ISTO  Experimental 
Facility  at  Cape  Canaveral;  further  analysis  awaits  more  accurate  information  on  the  system  point 
spread  function. 
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The  mathematical  formalism  underlying  these  restoration  algorithms  is  given  in  Appendix  A.  The 
algorithms  all  incorporate  regularization  to  counteract  the  destabilizing  effects  of  noise,  and 
require  an  estimate  of  the  object  support  (i.e.,  the  spatial  region  within  which  the  target  is  non¬ 
zero).  In  the  case  of  the  non-negative  reconstruction  algorithm,  this  estimate  is  refined  from  one 
iteration  to  the  next. 

(i)  Direct  recc  vistmction  from  the  regularized  pseudo  inverse. 

This  algorithm  is  widely  applicable.  It  can  be  applied  to  finely  sampled,  nearly  perfect  images  or 
to  highly  aberrated  shift-variant  images  in  which  pixellation  effects  are  clearly  evident.  The 
computati?n  is  direct  -  i.e.,  non-iterative  -  and  the  resulting  solution  can  exhibit  both  positive  and 
negative  values.  The  reconstruction  is  bounded  by  the  estimate  of  the  object  support  and  high 
resolution  can  be  achieved.  This  method  can  give  particularly  good  results  when  the  unresolved 
targets  of  interest  lie  within  a  limited  dynamic  range.  Weighting  can  be  incorporated  to  reflect  the 
level  of  confidence  in  the  data  and  varying  degrees  of  smoothness  can  be  enforced  on  the 
reconstruction. 

(a)  Unweighted  reconstruction 

This  practical  example  will  demonstrate  the  use  of  the  regularized  pseudoinverse  for  image 
reconstruction  with  the  identity  matrix  as  the  constr^t  operator  (see  Appendix  A).  Three 
small  incoherently-illuminated  apertures  were  imaged  onto  a  CCD  array  and  the  optical 
geometry  arranged  so  that  the  images  of  these  objects  all  lay  within  one  half  the  Rayleigh 
distance  of  each  other.  It  can  be  seen  from  Fig.  1  that  the  individual  images  are 
unresolved.  A  good  estimate  of  the  system  point  spread  function  is  essential  for  effective 
restoration,  and  I^g.  2  shows  the  result  of  avera^g  16  estimates  of  the  ps^  obtained  with 
a  point  light  source,  followed  by  filtering  in  the  Fourier  transform  plane  to  eliminate  out- 
of-band  noise.  The  ima^ng  aperture  was  constructed  with  a  pair  of  crossed  slits;  the  psf 
was  therefore  a  product  of  two  sine  functions.  Fig.  3  is  the  computed  reconstruction  of 
the  target  field.  The  objects  are  now  clearly  resolved.  Finally,  Fig.  4  shows  the  image 
obtained  with  a  large-aperture  lens  in  place  of  the  crossed  slits.  Contour  plots  of  Figs.  3 
and  4  show  that  the  reconstructed  images  are  in  their  correct  relative  locations. 
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(Sampling  Approximately  4x  Nyquist  Rate) 


Figure  2.  Filtered  Averaged  Point-Spread  Function 


Figures.  Reconstruction 


Figure  4.  High  resolution  Image  of  Object 


(b)  Weighted  reconstruction 


Weighting  can  be  introduced  into  the  calculations  in  order,  for  example,  to  discriminate 
against  bad  pixels  in  the  array,  or  to  adjust  the  relative  significance  in  the  computation  of 
bright  and  faint  features  in  the  target  area.  In  the  general  case,  the  weight  matrix  is 
derived  firom  a  covariance  matrix  which  expresses  information  on  the  correlations  between 
the  statistical  errors  in  the  image  data  values,  as  well  as  on  their  magnitude.  If  the  pixel 
weightings  are  uncorrelated,  the  weight  matrix  is  diagonal.  This  will  be  the  case  for  an 
image  in  which  fluctuations  ('shot  noise')  due  to  the  Poisson  statistics  of  the  detection 
process  are  the  dominant  noise  mechanism. 

(ii)  Reconstruction  of  subpixel  detail  fi’om  combinations  of  defocused  images 

In  addition  to  intrinsic  noise,  a  CCD  detector  degrades  the  image  information  in  two  ways.  First, 
there  is  necessarily  a  modest  amount  of  averaging  performed  over  the  active  area  of  the  pixel,  and 
second,  the  pixel  structure  imposes  a  spatial  sampling  rate.  When  dealing  with  CCD  data  sampled 
at  the  Nyquist  rate  of  the  imaging  system  the  degradation  caused  by  the  detector  is  not  great. 
Practically  speaking,  when  imaging  thermal  (incoherent)  radiation  one  requires  approximately  5 
samples  across  the  point  spread  function  core  to  achieve  the  Nyquist  rate.  When  the  optical 
design  results  in  coarser  sampling  (for  example,  when  the  point  spread  function  is  matched  to  the 
pbcel  size)  the  detector  itself  severely  limits  the  spatial  resolution  attainable.  Noise  fi'om 
background  sources  or  the  detector  elements  themselves  now  limits  the  useful  information  to  a 
relatively  few  pixels.  As  a  consequence,  although  estimates  can  be  made  of  the  centroid  of 
isolated  point  targets  to  subpixel  precision,  a  single  image  cannot  be  made  to  yield  detailed 
information  about  more  complex  object  fields. 

The  discussion  in  Section  3.1  of  the  image  restoration  problem  can  be  generalized  and  a  method 
devised  for  extracting  subpbcel  detail  fi’om  a  small  number  of  independent  images.  These  can  be 
generated  in  a  number  of  ways,  such  as  by  simultaneously  combining  defocused  and/or  laterally- 
shifted  images  of  the  same  cluster  of  targets.  One  could  also  use  time-series  data,  although  the 
evolution  of  the  image  of  a  moving  object  over  time  may  limit  the  resolution  obtainable. 
Alternatively  one  could  scan  an  image  with  a  1-D  array  so  that  high  sampling  rates  are  possible  in 
the  scan  direction,  and  relatively  low  rates  in  the  array  direction.  It  is  assumed  that  a  pixel  may 
contain  images  of  more  than  one  target,  which  poses  a  more  difficult  problem  than  that  of  locating 
the  position  of  a  single  target  with  sub-pixel  accuracy.  A  combination  of  correctly  chosen 
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defocused  images  has  been  found  to  produce  superior  results  to  a  combination  of  translated 
images  [12].  The  signal-to-noise  requirements  are  more  demanding  than  when  pixel  size  is 
assumed  to  represent  the  fundamental  limit  on  resolution. 

(iii)  Reconstruction  with  enforced  non-negativity 

Passive  thermal  imaging  systems  yield  non-negative  images,  and  this  constraint  may  be  important 
if  the  target  area  is  relatively  faint.  Greater  resolution,  suppression  of  ringing  and  larger  dynamic 
ranges  are  possible  when  a  non-negativity  constraint  is  imposed.  It  can  be  incorporated  by 
iterating  the  computation  a  limited  number  of  times;  pbcels  with  negative  content  can  be 
eliminated  after  each  iteration  or  discriminated  against  by  weighting.  In  this  manner  one  can 
adaptively  obtain  a  very  tight  target  support  constraint  with  little  a-priori  information.  Achieving 
tight  support  constraints  gives  the  greatest  degree  of  spectral  extrapolation,  and  therefore 
resolution.  Indeed,  for  small,  localizable  targets,  the  support  estimate  contains  most  of  the 
desired  targeting  information. 

Figs.  5-9  illustrate  the  performance  of  the  non-negative  algorithm  in  which  negative  pixels  are 
eliminated  at  each  iteration,  in  combination  with  the  method  for  recovering  sub-pixel  detail 
described  above.  Fig.  S  shows  the  object  field  superimposed  on  nine  pbcels  of  a  simulated  CCD 
device,  and  Figs.  6  and  7  their  images  on  the  central  seven  by  seven  array  at  varying  degrees  of 
defocus.  These  four  images  are  combined  in  sequence,  and  a  composite  matrix  formed  fi-om  the 
matrices  representing  the  corresponding  imaging  operators.  An  iterative  computation  using  the 
regularized  pseudoinverse  of  this  matrix,  with  non-negati\dty  enforced  at  each  step,  was  then 
carried  out.  The  result  for  the  'noiseless'  case  is  shown  in  Fig.  8.  Note  that  the  optimum  value  of 
the  regularization  parameter  is  non-zero  even  in  the  absence  of  added  noise;  computer  roundoff 
error  alone  is  sufficient  to  induce  instability  in  the  reconstruction.  It  can  be  seen  that  the 
individual  sources  have  been  correctly  recovered  in  number  and  position.  Fig.  9  shows  the  result 
obtained  when  Gaussian  noise  with  a  standard  deviation  of  5%  of  the  peak  pbcel  value  was  added 
to  the  data.  One  or  two  small  artifacts  have  appeared  and  there  is  some  blurring,  but  the 
individual  targets  are  clearly  identifiable. 


18 


3 

Z5 

2 

1 

OJ 

OJ  1  1.5  2  3 


■  1 

♦ 

!=■ - 

♦ 

1 - 

♦ 

_ ^ _ , 

♦ 

♦ 

♦ 

♦ 

_ 1 _ 

❖ 

1 

♦ 

1^ _ 1 

Figure  5.  Object  Field  on  Pixel  Structure 
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Figured.  Images  at  0  ana 
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Figure  7.  Images  at  1.6  and 2.4  Waves  Defocus 


Figure  8.  Reconstruction  with  OH  Added  Noise;  beta  =  5e-5 
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Ftgure  9.  Reconstruction  with  5H  Added  Noise;  beta  =  le-3 
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3.3  Computational  Aspects 


An  effective  image  restoration  algorithm  for  SDI  applications  must  not  only  exploit  all  available 
information  about  the  target,  it  must  also  be  well  matched  to  the  computational  hardware  and  be 
carefully  implemented  to  achieve  optimum  performance  in  near-real  time.  In  implementing  the 
algorithms  described  above,  two  difSculties  become  apparent  as  image  size  grows.  First,  the 
sheer  number  of  computations  can  be  very  large,  resulting  in  long  computing  times.  Second,  the 
matrix  sizes  in  two  dimensions  quickly  become  overwhelnungly  large.  For  example,  a  64  x  64 
image  reconstructed  into  64  x  64  elements  has  an  ima^g  operator  matrix  of  size  64^  x  64^,  or 
16.777  X  10®  elements.  For  single  precision  operation  (4  bytes/element),  this  results  in  67  Mbytes 
of  data  (134  Mbytes  with  double  precision).  Handling  such  large  data  sets  for  modest  images 
makes  the  high  performance  techniques  described  earlier  challenging  to  implement  even  for  earth- 
based  applications,  where  extensive  resources  are  available  and  time  may  not  be  a  factor.  Thus 
storage  requirements  and  computational  complexity  are  two  aspects  of  algorithm  implementation 
which  must  be  addressed  for  all  but  the  smallest  images. 

3.3.1  Reducing  Storage  Space 

There  are  essentially  three  approaches  to  reducing  space  demands.  The  first,  limited  to  shift 
invariant  systems,  is  to  use  an  iterative  procedure  based  on  the  fast  Fourier  Transform  (FFT). 
One  can  implement  the  regularized  pseudoinverse  reconstruction,  including  non-negativity,  in  this 
manner.  The  procedure  is  then  a  regularized  form  of  the  Gerchberg-Papoulis  algorithm  [13].  The 
data  and  image  sets  are  the  same  size,  and  only  FFTs  are  involved.  Hence  each  iteration  is  fast 
and  easily  implemented;  however,  convergence  is  slow,  and  1000  or  more  iterations  may  be 
required.  The  relative  performance  of  this  algorithm  has  not  been  evaluated. 

The  second  method  is  to  implement  the  regularized  pseudoinverse  exactly  by  partitioning  the 
matrix  equation  into  submatrices  and  using  a  sequential  calculation  based  on  smaller  matrices.  It 
has  been  shown  that  a  real  speed  gain  is  possible  for  the  pseudoinverse,  although  the  expression 
for  the  regularized  pseudoinverse  has  still  to  be  explored.  Smaller  matrices  are  manipulated,  and 
if  the  point  spread  function  is  shift  invariant,  some  of  the  submatrices  will  be  duplicated.  Non¬ 
negativity  probably  cannot  be  implemented  in  this  way. 

The  third  method  is  to  break  the  problem  into  smaller  blocks  and  deal  with  each  of  the  blocks 
independently  [14].  The  procedure  is  iterative,  and  potentially  faster  than  solving  the  original 
problem  as  a  unit.  However,  one  must  perform  additional  calculations  to  remove  edge  effects 
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from  neighboring  blocks.  (The  second  method  includes  all  cross-terms.)  The  ultimate  efficiency 
of  this  method  depends  on  the  number  of  iterations  required. 

3.3.2  Reducing  Computational  Complexity 

Considerable  effort  has  been  devoted,  in  collaboration  vdth  members  of  NOSC  staf^  to 
developing  techniques  for  alleviating  the  potential  computational  burden  associated  with  the 
algorithms  discussed  above  [IS,  16].  If  the  image  restoration  matrix  has  to  be  determined  for 
each  image,  the  number  of  operations  involved  in  computing  the  restoration  will  be  0(N^  for  an 
N  X  N  image.  For  even  modest  values  of  N  ,  this  number  will  be  very  large;  for  example,  if 
N  =  32,  the  number  of  operations  is  O(IO^)  . 

The  entries  in  the  ima^g  matrix  A  consist  of  discrete  values  of  the  system  point  spread 
function.  These  elements  are  determined  by  the  locations  of  the  samples  in  the  image  plane  and 
the  locations  of  the  points  at  which  the  corrected  values  are  to  be  determined  in  the  reconstruction 
plane.  The  point  spread  function  will  in  general  be  locally  shift-invariant,  but  change  gradually 
across  the  field  of  view,  and  it  will  be  assumed  that  the  instrument  is  fiiUy  characterized  by  a 
stored  set  of  calibration  data.  If  the  point  spread  function  is  also  time-varying,  recalibration  must 
be  performed  as  often  as  necessary  with  some  appropriate  source;  if  they  are  suitably  equipped, 
one  of  the  'eyes'  could  possibly  fulfill  this  function.  Thus,  for  a  given  target  area,  we  can  assume 
that  the  required  point  spread  function  is  already  available.  Suppose  that  image  size  and 
reconstruction  support  (the  region  within  which  the  target  objects  are  estimated  to  lie),  or  a 
typical  set  of  them,  are  known  in  advance.  Then  only  a  value  for  the  regularization  parameter  is 
required,  and  it  is  quite  likely  that  a  small  discrete  set  of  these  will  suffice.  (It  may  be  noted  here 
that  truncating  the  singular  values  of  A  at  the  appropriate  point  can  be  shown  to  be  equivalent  to 
incorporating  a  regularization  parameter  in  the  calculation  [17],  and  may  offer  some  advantages  in 
terms  of  parallel  processing  efficiency.)  Then  a  set  of  restoration  matrices  can  be  precomputed 
and  the  processing  of  the  image  reduced  to  a  single  vector-matrix  multiplication  involving  0(N^) 
operations.  This  can  still  be  a  daunting  task;  for  a  100  x  100  image,  0(10*)  operations  are 
required.  We  therefore  look  for  special  mathematical  features  in  the  structure  of  the  problem 
which  may  be  exploited  to  accelerate  the  computation  further. 

The  imaging  operator,  from  which  A  is  derived,  is  a  convolution  for  the  shift-invariant  system 
and  has  natural  cyclic  properties;  in  the  case  of  equal  numbers  of  samples  and  equal  sampling 
intervals  in  image  and  reconstruction  spaces,  the  matrix  A  is  Toeplitz.  Even  if  these  conditions 
are  not  met,  A  may  retain  much  Toeplitz-like  structure.  The  concept  of  displacement  rank  [18] 
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may  be  used  to  exploit  this  property;  the  singular  value  decomposition  of  A  then  provides  a 
means  for  constructing  approximations  to  the  image-restoration  matrix  (the  regularized  pseudo¬ 
inverse  of  A  )  in  the  form  of  relatively  low-order  sums  of  triangular  Toeplitz  matrix  products. 
Fast  multiplication  of  the  image  vector  by  the  image-restoration  matrix  then  becomes  possible 
[19].  Further  theoretical  developments  have  also  recently  appeared  in  the  open  literature  [20,  21]. 

For  the  two-dimensional  image,  the  imaging  matrix  A  can  be  written  as  a  matrix  consisting  of 
Toeplitz  sub-matrices  arranged  in  Toeplitz  order  (block-Toeplitz  with  Toeplitz  blocks).  Then  it 
can  be  shown  that  the  restoration  matrix  R  (the  regularized  pseudoinverse  of  A  )  is  block- 
persyirunetric  with  persymmetric  blocks.  Thus  it  retains  marked  cyclic  features.  On  the  basis  of 
numerical  experiments,  it  has  been  proposed  that  a  circulant  approximation  could  be  made  to 
R  ,  which  would  vastly  accelerate  the  final  computation,  since  FFT  techniques  would  now 
become  available.  Recent  work  [22]  in  this  field  provides  supporting  evidence  that  under  cert^ 
conditions  a  series  expansion  may  exist  for  R  in  which  the  first  term  is  a  circulant.  The  real 
significance  of  the  proposal  is,  of  course,  that  the  essential  information  needed  to  carry  out  the 
reconstruction  is  contained  in  a  single  row  or  colunm  of  R  .  A  further  possibility  now  suggests 
itself,  that  there  may  be  a  more  direct  way  to  obtain  this  information  than  via  the  QR  or  singular- 
value  decomposition  of  A .  In  the  most  recent  phase  of  this  program,  the  collaborative  research 
between  Spectron  and  NCCOSC  has  concentrated  on  this  topic. 

Displacement  rank  and  related  expansions  are  naturally  suited  to  algorithm  implementation  on 
parallel  processors,  and  the  latter  has  been  kept  constantly  in  mind  in  the  course  of  this  program. 
The  heavy  computational  demands  associated  with  image  reconstruction  Avill  be  most  effectively 
satisfied  by  a  judicious  combination  of  advanced  techniques  of  matrix  algebra  and  a  carefully- 
matched  parallel-processing  architecture. 
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4.0  WIDER  APPLICATIONS 


The  key  to  the  applicability,  for  any  ^ven  task,  of  the  algorithms  discussed  in  this  document  is 
that  the  "image"  and  the  "object"  are  related  by  a  linear  transformation.  Then,  for  a  sampled  data 
system,  the  discrete  representation  of  the  image  (the  sensor  data)  and  the  discrete  approximation 
to  the  object  will  be  related  by  a  matrix.  The  matrix  incorporates  all  the  evolutionary  processes 
undergone  by  the  information-carrying  field  during  propagation  from  the  object  and  collection, 
detection  and  digitization  by  the  sensor  system.  In  the  case  of  a  conventional  optical  image,  the 
field  is  electromagnetic,  and  the  transformation  consists  primarily  of  Fourier  transforms;  the 
sensed  data  may  be  further  modified,  of  course,  by  the  characteristics  of  the  detector  and  the 
associated  circuitry.  Nonlinear  effects,  such  as  Poisson  noise  in  the  optical  case,  can  also  be 
accommodated  provided  they  are  not  so  great  that  the  assumption  of  a  linear  object-image 
relationship  ceases  to  be  a  reasonable  approximation. 

Potential  applications  exist  in  many  areas  of  remote  sensing,  including  tomography,  synthetic 
aperture  radar,  ima^g  radar,  sonar,  seismic  prospecting,  spectroscopy  and  system  identification 
and  control.  In  the  imaging  radar  case,  for  example,  tomographic  reconstruction  of  the  scattering 
function  via  the  squared  modulus  of  the  transmitted  signal's  ambiguity  function  (the  analog  of  the 
optical  point  spread  function)  would  yield  enhanced  resolution  [23].  In  the  case  of  synthetic 
aperture  radar,  the  point  spread  function  is  separable  and  particularly  efficient  reconstructions  of 
the  scattering  distribution  are  possible.  It  may  also  be  noted  that  missing  data  segments,  which 
are  a  common  feature  of  inverse  scattering  problems,  can  be  recovered  with  the  techniques 
described  here. 
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5.0  CONCLUSIONS  AND  RECOMMENDATIONS 


As  a  result  of  the  major  geopolitical  events  which  have  recently  occurred,  the  nature  of  the 
ballistic  missile  attack  threat  to  the  United  States  has  undergone  significant  change.  Many 
additional  areas  of  the  globe  have  become  potential  launch  sites,  and  the  defense  of  regions 
outside  the  United  States  is  now  declared  public  policy.  The  surveillance  function  is  a 
fundamental  component  of  defensive  strategy,  and  a  clear  implication  of  these  developments  is 
that  a  comprehensive  network  of  space-based  observation  platforms  will  be  needed  to  meet  the 
SDI  requirements  for  detection,  discrimination  and  subsequent  tracking.  In  fulfilling  the 
surveillance  mission,  the  performance  of  the  on-board  optical  sensors  and  the  quality  of  the 
images  they  create  will  be  of  paramount  importance. 

Given  adequate  illumination,  the  performance  of  an  optical  system  is  limited  ultimately  by  the 
physical  laws  governing  the  phenomenon  of  dif&action.  Further  limitations  can  arise  from  design 
errors,  faulty  construction,  the  stresses  of  launch  and  deployment  or  the  effects  of  extended 
operation  in  the  space  environment.  This  report  has  been  concerned  with  the  exploitation  of 
modem  mathematical  techniques,  in  combination  with  state-of-the-art  computational  technology, 
to  compensate  for  such  defects  and  deficiencies.  The  feasibility  of  recovering  target  information 
apparently  irretrievably  lost  in  the  process  of  forming  and  capturing  the  image  has  been 
demonstrated  with  numerous  simulations  as  well  as  experimentally-derived  data.  As  described 
above,  algorithms  have  been  designed  and  developed  for  aberration  correction,  to  achieve 
resolution  beyond  the  diffiaction  limit  and  for  other  specific  image-enhancing  purposes  which 
require  only  the  provision  of  appropriate  computing  resources.  Powerfiil  techniques  of  linear 
algebra  exist,  including,  for  example,  displaced-difference  expansions  and  circulant 
approximations,  which  can  be  exploited  to  reduce  the  associated  computational  burdens.  Further 
research  is  needed  in  this  area  into  the  ranges  of  validity  of  the  approximations  and  their  most 
efiBcient  implementation. 

The  possible  problems  which  may  arise  fi-om  the  stresses  of  launch  and  deployment,  fi'om  errors 
of  manufacture  which  may  not  be  discovered  until  the  system  is  in  orbit  and  operational,  and  fi'om 
aging  of  the  optical  and  electronic  components,  imply  that  provision  for  recalibration  of  the 
optical  subsystem  must  be  made.  This  can  most  easily  be  accomplished  using  a  cooperative  target 
equipped  with  a  light  source,  which  could  be  another  eye,  if  these  are  active.  Thus 
communications  channels  between  the  orbiting  eyes  will  be  needed.  In  addition,  computation  of 
the  modified  reconstruction  matrix  and  subsequent  updating  of  the  on-board  data  would  be  most 
efficiently  performed  on  a  large  ground-based  computer,  requiring  a  ground-orbit  communication 
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link  for  downloading  of  the  calibration  data  and  uploading  of  the  new  matrix  elements.  State-of- 
the-art  computer  technology  would  be  available  for  the  near  real-time  correction  of  system 
deficiencies  on  a  continuing  basis,  thereby  considerably  extending  the  useful  lifetime  of  the 
orbiting  sensor  network. 

It  is  therefore  recommended  that,  in  addition  to  continuation  of  the  research  into  algorithm 
structure  and  implementations,  provision  be  made  for  both  adequate  computing  resources  and 
appropriate  sensor-to-sensor  and  ground-to-sensor  conununication  links  in  the  basic  design  of  the 
surveillance  sensor  system,  thus  making  possible  both  the  correction  of  image  deficiencies  and  the 
efficient  updating,  as  necessary,  of  the  computational  data  base.  The  benefits  will  include 
enhanced  detection,  discrimination  and  tracking  capabilities  and  a  longer-lived,  more  cost- 
effective  defensive  screen. 
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APPENDKA 


A-1 


REGULABIZED  IMAGE  RECONSTRUCTION 


The  task  of  image  reconstruction  belongs  to  the  more  general  class  of  ill-posed  inverse  problems, 
which  arise  in  many  fields  besides  optics  -  for  example  in  radar,  geophysics,  tomography  or  the 
extrapolation  of  band-limited  signals.  Suppose  we  are  given  the  operator  equation 

Ax  =  y 

where  A:  X->Y .  The  problem  of  finding  the  solution  x,  given  y,  is  said  to  be  well-posed  if 
for  each  y  there  exists  a  unique  solution  which  depends  continuously  on  the  data;  i.e.,  is  stable 
under  small  changes  in  the  data  [G].  Otherwise,  the  problem  is  ill-posed.  It  should  be 
emphasized  that  ill-posedness  is  inherent  in  the  physical  nature  of  the  measurement  process,  not 
the  mathematical  modeling  of  it. 

Consider  a  typical  image  reconstruction  problem.  We  suppose  some  object  f ,  an  element  of 
Hilbert  space  F ,  is  imaged  by  a  linear  operator  A  into  an  element  of  another  Hilbert  space  G , 
and  that  this  image  is  perturbed  by  an  unknown  additive  noise  fimction  which  may  arise  in  dther 
the  imaging  or  the  observation  process.  The  effect  of  the  noise  in  practice  will  generally  be  to 
displace  the  image  g  out  of  the  range  of  i4  -  there  will  no  longer  be  an  element  of  F 
corresponding  to  g  and  the  equation  Af  =  g  will  have  no  solution.  Clearly  the  problem  is  ill- 
posed. 

This  difficulty  of  a  lack  of  continuous  dependence  can  be  overcome  by  one  of  the  various 
techniques  of  "regularization."  Essentially,  the  ill-posed  problem  is  replaced  by  a  related  well- 
posed  one^  chosen  to  be  physically  meaningful  and  to  possess  the  necessary  properties  of 
convergence  and  stability.  Thus  we  may  change  the  concept  of  a  solution,  or  the  Hilbert  spaces 
or  their  topologies,  or  the  operator  itself  The  technique  we  shall  use  belongs  to  the  last  category. 

We  impose  physically  reasonable  constraints  on  the  permitted  solutions.  If  £  is  a  measure  in  the 
norm  sense  of  the  noise  in  the  image  (norm  in  the  appropriate  Hilbert  space  is  denoted  by 
II  •  1^  )  and  if  C  is  some  constraint  operator  with  E  a  known  bound,  we  shall  require  that  all 
possible  reconstructions  f  satisfy 

\g-^f'\o  ^  ^  and  ||C/'||, 
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C  may  be  used,  for  example,  to  impose  smoothness  on  the  reconstruction  or  to  weight  the 
reconstruction  support.  If  C  is  the  identity  operator,  £  is  a  bound  on  the  norm  of  the 
reconstruction.  We  combine  the  constraints  quadratically  and  minimize  the  functional 


where  fi  =  i'fl^  .  Note  that  smaller  values  of  p  are  equivalent  to  demanding  greater  fidelity 
between  the  reconstruction  and  the  data;  greater  values  place  more  emphasis  on  the  property  of 
the  reconstruction  associated  with  C.  In  the  present  discussion  we  shall  take  C  =  I.  The 

minimizer  can  be  expressed  in  either  of  the  forms 

f,  =  {A>A  +  fil)-'  A-g 
or 

/,  =  A>(AA»  +  fil)-'  g 

where  A*  is  the  operator  adjoint  to  A.  The  inverses  of  the  bracketed  operators  will  always 
exist,  since  the  eigenvalues  of  the  symmetric  operators  A  *A  and  AA*  are  non-negative. 

The  operator,  the  image  and  the  reconstruction  will  consist  in  practice  of  finite  arrays.  We  shall 
assume  that  the  imaging  operator  has  been  recast  necessary)  as  a  two-dimensional  matrix  and 
that  /  and  g  are  vectors.  The  finite-dimensional  approximations  to  the  Hilbert  spaces  F  and  G 
are  ^  and  ,  N~  and  A/-dimensional  spaces  of  real  or  complex  numbers.  For  matrices  with 
complex  entries.  A*  becomes  A^  and  the  solution  is 

4  =  Ar  ,  (I) 


where  R  can  take  either  of  the  forms 
(a* A  +  fil,)-'  A-g 
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or 

^"(^4*  +  g  . 

Either  of  these  matrix  expressions  is  referred  to  as  the  regularized  pseudoinverse  of  ^4.  To 
facilitate  analysis  and  gtun  insight  into  the  nature  of  the  ill-posed  problem  we  make  use  of  the 
method  of  singular  value  decomposition  (SVD).  (There  is  also  a  significant  gain  in  accuracy  in 
using  the  SVD  rather  than  in  forming  A^A  directly.)  We  write  for  the  SVD  of  A 

A  =  imV“  (2) 

where 

U»u  =  V”V  =  W“  =  1^ 
and 

Z  =  <*urg(oi....crJ,  a,  ^  0. 

U  consists  of  the  N  othonormalized  eigenvectors  ii,  associated  with  the  N  largest  eigenvalues 
of  AA^^  and  V  consists  of  the  N  orthonormalized  eigenvectors  v,  of  A^A.  The  Oj ,  the 
non-negative  square  roots  of  the  eigenvalues  of  A^A ,  are  the  singular  values  of  A  . 

Then,  fi'om  equation  (1), 

If  ^  =  0  ,  we  have 

f,=A*g 

where  A*  denotes  the  Moore-Penrose  generalized  inverse  (familiarly  called  the  pseudoinverse) 
of  A  .  Thus  the  matrix  operator  in  equation  (1)  can  accurately  be  described  as  a  regularized 
pseudoinverse  of  A. 

Tlw  regularizing  role  of  p  can  clearly  be  seen  in  equation  (3);  for  zero  p ,  the  elements  of  the 
diagonal  matrix  are  simply  the  inverses  of  the  non-zero  singular  values  of  A  and  can  grow 
without  bound  as  these  tend  toward  zero. 
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Weighting  can  be  introduced  into  both  of  the  spaces  F  and  G  through  the  definition  of  the  inner 
product,  and  can  be  used  to  predispose  the  shape  of  the  reconstruction  to  some  prior  expectation, 
or  to  reflect  statistical  information  about  the  noise.  A  modified  form  of 
equation  (1),  incorporating  the  weighting  matrices,  now  governs  the  reconstruction,  which  with 
the  aid  of  Choleslqr  factorization  can  again  be  analyzed  by  SVD  [D]. 
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ILLUSTRATIVE  IMAGE-ENHANCEMENT  PROGRAMS 


A  suite  of  operational  programs  based  on  selected  image-enhancement  algorithms  which  were 
developed  under  this  contract  can  be  found  on  the  floppy  disk  included  herewith.  The  programs 
are  designed  to  run  in  the  MATLAB  environment  on  a  Sun  workstation.  Their  individual 
functions  are  sununarized  below.  The  reconstruction  procedures  are  based  on  the  assumption  that 
the  point-spread  function  (psf)  is  two-dimensional  and  non-separable. 

A  demonstration  is  available  by  typing  »startup2d.  This  routine  generates  and  saves  the  psf  and 
the  singular-value  decomposition  of  the  psf  for  the  demonstration  procedure.  For  convenience, 
(sinc(x)*sinc(y)).^2  was  used  for  the  psf  This  function  is,  in  fact,  separable,  although  the 
computation  is  performed  as  if  it  were  non-separable.  Note  that  'infinities'  are  generated  in 
attempting  to  calculate  sin(0)/(0) ;  these  result  in  standard  MATLAB  warning  messages,  but  have 
no  effect  on  the  reconstruction. 

demo2d_4  is  called  at  the  end  of  the  startup2d  routine;  after  the  first  setup  run,  it  can  be  called 
directly.  The  program  then  randomly  constructs  an  object,  forms  its  image,  adds  the  specified 
amount  of  noise  and  calculates  a  value  for  the  regularization  parameter  using  the  technique  of 
weighted  cross-validation.  Finally  the  reconstruction  is  computed,  using  two  different  methods 
for  estimating  the  object  support.  The  results  are  displayed  in  two  final  graphics  screens  (note  the 
<PAUSE>  between  them). 

Reasonably  accurate  estimation  of  the  object  support  is  an  important  requirement  if  good 
reconstructions  of  an  object  are  to  be  achieved.  In  this  demonstration,  the  first  method  for 
support  estimation  uses  a  reconstruction  space  which  is  coextensive  with  the  space  over  which  the 
object  was  generated.  This  constrrunt  is  very  weak  -  it  would  be  somewhat  better,  for  example,  to 
use  a  smaller  bound  estimated  fi'om  the  image  width  and  the  psf  width  -  and  the  resulting 
reconstruction  is  often  not  a  great  improvement  over  the  original  image.  However,  the  second 
algorithm  (niirec_w2)  generally  performs  extremely  weU.  In  estimating  the  reconstruction 
support,  this  algorithm  uses  the  fiict  that  the  object  should  be  non-negative.  It  calculates  the 
reconstruction  in  10  iterations  (the  succession  of  2-d  reconstructions  mapped  column-by-column 
to  vectors  is  displayed  during  this  process).  The  assumption  is  made  after  each  step  that  any 
negative  points  imist  be  outside  the  object  support,  and  a  new  reconstruction  is  generated  with  the 
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modified  support.  This  procedure  is  most  powerful  in  the  2-d  nonseparable  case,  and  cannot  be 
efiectively  applied  to  2d-separable  problems  where  the  separability  is  exploited  in  the 
computations,  since  changing  one  pixel  in  one  of  the  dimensions  will  in  that  case  affect  an  entire 
row  of  the  image. 


The  figures  at  the  end  of  this  Appendix  are  an  example  of  the  output  of  the  program  demo2d_4. 


Program  Function  Summary 


2-D  Mapping  Algorithms 

mat2vec  Takes  a  matrix  (such  as  an  image)  and  maps  it  to  a  vector. 

vec2mat  Takes  a  vector  and  maps  it  to  a  matrix. 

mul  A  convenient  procedure  for  multiplying  A  and  f,  where  A  is  the  psf  matrix 

and  f  is  in  the  form  of  a  2-d  image,  mul  maps  f  to  a  vector,  performs  the 
multiplication,  and  remaps  A*f  back  into  a  2-d  image. 


Point-Spread  Function  Definition 

mkps£2d_sinc  Forms  sinc^2psf  matrix.  This  psf  corresponds  to  a  square  unapodized 
aperture  (a  slit). 

mkpsQd_sincl2  :  Forms  sinc^2  psf  matrix  with  the  sampling  rate  in  image  space  equal  to  half 
that  of  object  space.  This  psf  corresponds  to  a  square  unapodized  aperture 
(asUt). 


Regularized  Matrix  Inversion 

oper  Forms  inv(A'A  +  beta*I)A'  using  MATLAB  inv  function. 

operpen  Forms  the  pseudo-inverse  of  A  from  the  svd  of  A. 

opersvd  Forms  regularized  pseudo-inverse  from  svd  of  A . 

operwf  Forms  inverse  of  inv(A'A+betaW^A' ,  where  W  is  specified  by  the  user, 

and  weights  the  reconstruction  norm. 
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recld 


Calculates  the  regularized  pseudo-inverse  Q  and  the  reconstruction.  Allows 
for  the  specification  of  a  convex  support  constraint. 


recldjg  Calculates  a  reconstruction  using  Gaussian  elimination.  Allows  for  the 

specification  of  a  convex  support  constraint. 

recld_s  Calculates  a  reconstruction  using  svd  of  psf  Allows  for  the  specification  of  a 

convex  support  constraint. 

nnrec _gl  Calculates  non-negative  reconstruction  using  Gaussian  elimination.  Effectively 

derives  support  constraint  in  the  process.  The  support  decreases  firom  iteration 
to  iteration,  speeding  the  calculation,  but  there  is  a  risk  of  eliminating  pbcels 
within  the  true  support. 

nnrec  jg2  Essentially  the  same  as  nnrec _gl,  except  that  the  regularization  parameter 

decreases  with  each  iteration. 

nnrec_wl  Calculates  a  non-negative  reconstruction  uung  Gaussian  elimination.  Uses 

inner  product  weighting  to  eff^vely  derive  a  support  constraint.  Differs  firom 
nnrec ^1  in  that  computational  time  per  iteration  is  constant,  but  the  support 
can  evolve  over  time. 

nnrec  w2  nnrec  wl  limited  to  10  iterations. 


xvalidl  Calculates  the  regularization  parameter  for  a  1-d  image. 

xvalid2  Calculates  the  regularization  parameter  for  a  block  of  1-d  images. 


xvalid2d_l  Calculates  the  cross-validation  fiinction  over  a  specified  range  of  the 
regularization  parameter  for  a  ^ven  image. 
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xvldj 

xvld_s 


Carries  out  the  cross-validation  Wahba  sum  using  the  svd  of  A. 

Carries  out  a  truncated  sum  similar  to  xvldj,  but  generally  shows  increased 
variance  in  the  cross-validation  curve  for  small  values  of  the  regularization 
parameter,  for  1-d  images  (caused  by  insufBcient  data  points). 
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ABSTRACT 

Image  restoration  procedures  are  commonly  unstable  in  the  presence  of  noise,  and  some  technique 
for  restoring  stability  becomes  essential.  The  methods  of  regularization  theory  are  particularly 
appropriate  for  this  purpose.  A  specific  type  of  regularized  solution  is  introduced  in  the  general 
context  of  image  reconstruction.  A  super-resolution  problem  is  then  considered  from  the  point  of 
view  of  the  computational  tasks  involved,  with  particular  reference  to  the  estimation  of  certain  key 
parameters  and  to  implementations  which  increase  the  efficiency  of  the  calculations.  Parameter 
estimation  is  performed  by  weighted  cross-validation.  The  improvement  in  efficiency  is  achieved 
through  the  exploitation  of  symmetries  or  cyclic  properties  inherent  in  the  reconstruction  operator. 
The  concept  of  displacement  rank  is  introduced  and  estimates  made  of  the  computational  burden 
associated  with  various  classes  of  regularized  reconstruction  matrices. 

INTRODUCTION 

We  shall  be  concerned  with  the  construction  and  computation  of  solutions  to  the  problem  of  image 
reconstruction  from  noisy,  incomplete  or  otherwise  degraded  data,  when  the  imaging  operator  is 
known.  The  strongly  smoothing  operators  which  are  typical  of  imaging  or  measurement  systems 
are  associated  in  the  reconstruction  process  with  an  extreme  sensitivity  to  noise  or  distortion  in  the 
data;  small  changes  induced  in  the  image  can  result  in  grossly  differing  reconstructions.  The 
immediate  cause  of  instability  in  the  computational  process  is  the  ill-conditioning  of  the  matrix 
representing  the  imaging  operator. 

The  task  of  image  reconstruction  in  fact  belongs  to  the  more  general  class  of  ill-posed  inverse 
problems,  which  arise  in  many  fields  besides  optics  -  for  example  in  radar,  geophysics,  tomography 
or  the  extrapolation  of  band-limited  signals.  Suppose  we  are  given  the  operator  equation 

Ax  =  y 

where  A  :  X  Y.  The  problem  of  finding  the  solution  x,  given  y,  is  said  to  be  well-posed  if  for 
each  y  there  exists  a  unique  solution  which  depends  continuously  on  the  data;  i.e.,  is  stable  under 
small  changes  in  the  data^  Otherwise,  the  problem  is  ill-posed.  It  should  be  emphasized  that  ill- 
posedness  is  inherent  in  the  physical  nature  of  the  measurement  process,  not  in  the  mathematical 
modeling  of  it. 
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Consider  a  typical  image  reconstruction  problem.  We  suppose  some  object  f  ,  an  element  of  a 
Hilbert  space  F  ,  is  imaged  by  a  linear  operator  A  into  an  element  of  another  Hilbert  space  G  , 
and  that  this  image  is  perturbed  by  an  unknown  additive  noi.se  function  which  may  arise  in  either 
the  imaging  or  the  observation  process.  The  effect  of  the  noise  in  practice  will  generally  be  to 
displace  the  image  g  out  of  the  range  of  A  -  there  will  no  longer  be  an  element  of  F  corresponding 
to  g  and  the  equation  Af  =  g  will  have  no  solution.  Clearly  the  problem  is  ill-posed. 

This  difficulty  of  a  lack  of  continuous  dependence  can  be  overcome  by  one  of  the  various 
techniques  of  "regularization".  Essentially,  the  ill-posed  problem  is  replaced  by  a  related  well- 
posed  one,  chosen  to  be  physically  meaningful  and  to  possess  the  necessary  properties  of 
convergence  and  stability.  Thus  we  may  change  the  concept  of  a  solution,  or  the  Hilbert  spaces  or 
their  topologies,  or  the  operator  itself.  The  technique  we  shall  use  belongs  to  the  last  category. 


THE  REGULARIZED  RECONSTRUCTION 


We  impose  physically  reasonable  constraints  on  the  permitted  solutions.  If  £  is  a  measure  in  the 
norm  sense  of  the  noise  in  the  image  (norm  in  the  appropriate  Hilbert  space  is  denoted  by  1 1  •  1 1  ) 
and  if  C  is  some  constraint  operator  with  E  a  known  bound,  we  shall  require  that  all  possible 
reconstructions  /'  satisfy^ 

Mg-A/’Mc^e  and  MC/'llp^E 

C  may  be  used,  for  example,  to  impose  smoothness  on  the  reconstruction  or  to  weight  the 
reconstruction  support.  If  C  is  the  identity  operator,  E  is  a  bound  on  the  norm  of  the  reconstruction. 
We  combine  the  constraints  quadratically  and  minimize  the  functional 

Mg-A/'Mc  +  3  Me/’ I  Ip 

where  3  =  e^/E^.  Note  that  smaller  values  of  3  are  equivalent  to  demanding  greater  fidelity 
between  the  reconstruction  and  the  data;  greater  values  place  more  emphasis  on  the  property  of  the 
reconstruction  associated  with  C.  In  the  present  discussion  we  shall  take  C  =  I .  The  minimizer /p 
can  then  be  expressed  in  either  of  the  forms 

/p  =  (A*A  +  3l)-lA*g 

or  /p  =  A"  (AA*  +  pi)-Ig 

where  A*  is  the  operator  adjoint  to  A.  The  inverse  of  the  bracketed  operators  will  always  exist, 
since  the  eigenvalues  of  the  symmetric  operators  A*A  and  AA*  are  non-negative. 
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Weighting  can  be  introduced  into  both  of  the  spaces  F  and  G  through  the  definition  of  the  inner 
product,  and  can  be  used  to  predispose  the  shape  of  the  reconstruction  to  some  prior  expectation,  or 
to  reflect  statistical  information  about  the  noise.  A  modified  form  of  equation  (1),  incorporating  the 
weighting  matrices,  now  governs  the  reconstruction,  which  with  the  aid  of  Cholesky  factorization 
can  again  be  analyzed  by  SVD^. 

PARAMETER  ESTIMATION  USING  WEIGHTED  CROSS-VALIDATION: 

To  perform  the  computations  in  equation  (1)  or  (3),  we  need  estimates  of  the  regularization 
parameter  P  and  of  the  region  in  object  space  (the  object  support)  into  which  the  reconstruction 
will  be  made.  In  arriving  at  these  estimates,  the  fullest  use  should  be  made  of  available  a  priori 
information  about  the  object  and  the  data-gathering  process.  Weighted  cross-validation^  provides  a 
method  for  estimating  a  near-optimum  value  of  the  regularization  parameter  from  the  data  alone. 
We  demonstrate  its  use  for  this  purpose  and  also  for  the  estimation  of  the  reconstruction  support  in 
the  super-resolution  problem. 

ESTIMATING  THE  REGULARIZATION  PARAMETER 

It  is  possible  to  obtain  a  good  estimate  of  the  optimum  regularization  parameter  from  the  data  alone 
through  the  use  of  weighted  cross-validation.  The  great  merit  of  this  technique  is  that  it  requires  no 
knowledge  of  the  object,  and  places  very  loose  constraints  on  the  image  and  noise,  namely  that  the 
noise  should  be  white  and  that  the  image  be  "very  smooth"  -  this  notion  is  defined  precisely  in 
Reference  5. 

The  idea  behind  weighted  cross-validation  is  as  follows.  Suppose  /  p  is  the  minimizer  of 

N  2 

5^  {(Ah)-gj}2  +  Pllhllp 

j=l 

j^k 

where  the  k^**  data  point  has  been  omitted.  We  can  now  predict  a  value,  (A/ p  ,  for  this  k***  point 
and,  for  a  given  value  of  P,  form  the  weighted  prediction  error  V(p): 

N 

V(P)  -  I  <(A/„j,)k-gk>2w„(p) 
k=l 

where  the  wj^CP)  are  weighting  functions  which  represent  the  relative  significance  of  the  g^  in  the 
prediction.  They  take  the  form 

Wk(p)=  r  (AAH  +  pi),-l  p 
L  trace  {(AA^  +  Pl)'* }  J 

It  can  be  shown  (Wahba)  that  the  P  which  minimizes  the  expression  above  for  V(P)  is  an  estimate 
of  the  minimizer  of  the  true  prediction  error: 
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MATRIX  OPERATORS  AND  THE  SVD 


The  operator,  the  image  and  the  reconstruction  will  consist  in  practice  of  finite  arrays.  We  shall 
assume  that  the  imaging  operator  has  been  recast  (if  necessary)  as  a  two-dimensional  matrix  and 
that  f  and  g  are  vectors.  The  finite-dimensional  approximations  to  the  Hilbert  spaces  F  and  G  are 
and  ,  N-  and  M-dimensional  spaces  of  real  or  complex  numbers.  For  matrices  with 
complex  entries.  A*  becomes  A^  and  the  solution  is 

/p  -  (AHA  +  PIn)->  AHg  , 

or  /p  -  AH(AAH-fpI„)-lg 

To  facilitate  analysis  and  gain  insight  into  the  nature  of  the  ill-posed  problem  we  make  use  of  the 
method  of  singular  value  decomposition  (SVD)^.  (There  is  also  a  significant  gain  in  accuracy  in 
using  the  SVD  rather  than  in  forming  A^A  directly.)  We  write  for  the  SVD  of  A 

A  =  ULVH  (2) 

where 

UHu  =  vHV  =  VVH  =  In 
and 

2  =  diag  (Oi . . .  On),  Oj  >  0. 

U  consists  of  the  N  orthormalized  eigenvectors  Uj  associated  with  the  N  largest  eigenvalues  of  AA^, 
and  V  consists  of  the  N  orthonormalized  eigenvectors  vj  of  A^A.  The  Oj,  the  non-negative  square 
roots  of  the  eigenvalues  of  A^A,  are  the  singular  values  of  A. 

Then,  from  equation  (1), 

/p  =  V  diag  (... ,  Oi/(Oi2  +  p) , ...)  uHg  (3) 

If  P  =  0 ,  we  have 

/p  =  A+g 

where  A'*'  denotes  the  Moore-Penrose  generalized  inverse  (familiarly  called  the  pseudo-inverse)  of 
A.  Thus  the  operator  in  equation  (1)  can  accurately  be  described  as  a  regularized  pseudo-inverse  of 
A. 

The  regularizing  role  of  P  can  clearly  be  seen  in  equation  (3);  for  zero  P,  the  elements  of  the 
diagonal  matrix  are  simply  the  inverses  of  the  non-zero  singular  values  of  A  and  can  grow  without 
bound  as  these  tend  toward  zero. 
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where  g  represents  the  ideal  noiseless  image 
The  estimator  VO)  can  be  expressed  in  the  form 

N  2 

I  {(AAH  +  pi)-lg}^. 

VO)  =  ^ _ 

[  trace{(AAH  +  pi)'*}]^ 

Thus  to  calculate  V(P)  it  is  necessary  to  invert  AA^  +  pi  .  In  the  general  case,  V(P)  can  be 
rewritten  in  terms  of  the  SVD  of  A: 

A  =  UZV** ,  L  =  diag  ( cTj , . . . ,  ). 

We  find 


VO) 


k=l  ^Ok^+p  n.2+R^  ^ 


Ok^+P 


where 


g  =  UHG 

Having  carried  out  the  SVD  of  A,  this  form  has  the  advantage  that  recalculating  V(p)  for  a  new 
value  of  P  involves  little  more  than  4N  multiplications. 

A  further  simplification  can  be  made  when  A  is  a  circulant  matrix.  In  this  case  Af  is  a  circular 
convolution,  and  the  calculation  can  be  carried  out  using  a  discrete  Fourier  transform  (DFT),  or, 
more  efficientiy,  a  fast  Fourier  transform  (FFT).  U  in  this  case  is  the  DFT  matrix  and  ^  is  the 
DFT  of  the  data. 

Although  the  problem  of  finding  the  regularization  parameter  would  now  be  highly  tractable,  this 
approach  is  not  usually  applicable  to  superresolution  for  two  reasons.  The  first  reason  is  that  in 
general  the  point-spread  function  in  an  imaging  problem  is  neither  a  periodic  function  nor  of  finite 
extent,  so  that  A  cannot  be  a  circulant.  However,  A  can  often  be  approximated  as  a  circulant  The 
second  consideration  is  relevant  in  spectral  extrapolation:  a  circulant  approximation  to  A  enforces 
the  same  support  bounds  on  the  reconstruction  as  on  the  data  g  ,  so  that  only  residual  extrapolation 
is  obtained  outside  the  system  bandpass.  It  might  be  possible  in  some  cases  to  alleviate  this 
problem  by  estimating  the  regularization  parameter  with  a  circulant  approximation  to  A  and  then 
carrying  out  the  actual  reconstruction  with  the  correct  A. 
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ESTIMATING  THE  OBJECT  SUPPORT 


A  number  of  approaches  are  possible  to  the  problem  of  estimating  the  object  support,  which  enters 
the  calculation  through  the  matrix  A.  Given  the  point-spread  function  of  an  imaging  system,  and 
possibly  prior  knowledge  of  the  types  of  objects  involved,  one  can  sometimes  form  a  good  estimate 
of  the  outer  boundary  of  the  object  support  simply  from  the  width  of  the  image.  An  alternative 
strategy  can  be  based  on  the  fact  that,  for  values  of  ^  around  the  optimum,  the  error  estimator  V(3) 
is  found  to  increase  rapidly  as  the  object  support  bound  is  reduced  below  the  correct  value.  (The 
analytical  behavior  of  V(P)  under  these  circumstances  is  still  under  investigation.)  This  suggests 
that  we  should  choose  the  support  bound  associated  with  the  minimum  of  the  surface  V(P,  support 
bound). 

An  image  with  5%  Gaussian  additive  noise  was  formed  of  the  object  shown  dashed  in  Figure  1. 
The  object  support  was  17  units  wide.  A  sequence  of  cross-validation  curves  was  calculated  over 
the  range  10'^<  P<  10^  for  reconstruction  supports  of  1  to  25  pixels.  The  resulting  surface  is 
shown  in  Figure  2.  The  minimum  occurs  at  the  point  corresponding  to  P  =  0.1,  support  =  17,  A  set 
of  values  of  V(P)  was  then  calculated  for  P  =  0.1  and  for  a  series  of  reconstruction  supports 
consisting  of  the  correct  support  bound  of  19,  but  differing  by  the  deletion  of  successive  internal 
single  points.  Thus  the  reconstruction  was  forced  to  zero  on  one  of  the  17  pixels  in  each  case. 
Figure  3  shows  the  resulting  values  of  V(P)  plotted  against  the  corresponding  deleted  pixels.  The 
positional  dependence  approximately  conforms  to  the  correct  object  support 

DISPLACEMENT  RANK  AND  FAST  MATRIX- VECTOR  MULTIPLICATION 

Inversion  of  an  NxN  matrix  takes  in  general  OfN^)  operations.  If  the  matrix  is  Toeplitz  the  number 
of  operations  needed  falls  to  0(N2).  If  the  imaging  operator  is  shift-invariant  and  the  numbers  of 
samples  and  sampling  rates  (for  unit  magnification)  are  the  same  in  image  and  reconstruction 
spaces,  then  the  imaging  matrix  A  is  Toeplitz.  However,  it  is  often  the  case  that  A  is  not  Toeplitz, 
but  nevertheless  retains  much  Toeplitz-like  structure.  This  is  true,  for  example,  if  we  assume  that 
the  point-spread  function  is  shift-invariant,  but  that  A  is  rectangular,  or  that  the  sampling  rates  in 
image  and  reconstruction  spaces  are  unequal,  or  that  columns  are  deleted  from  A  to  restrict  the 
reconstruction  support  in  some  appropriate  way.  In  all  of  these  cases,  A  retains  some  regularity  in 
structure  which  one  would  expect  could  be  used  to  make  computations  involving  A  more  efficient. 

The  concept  of  displacement  rank  was  introduced  to  quantify  these  ideas^.  The  SVD  then  provides 
a  mechanism  for  constructing  approximations  to  A  in  a  form  in  which  its  special  properties  can  be 
exploited^.  We  define  first  the  NxN  downshift  matrix  Z: 


0  0  ...  o' 

r 

Zx  = 

1  0  ...  0 

X2 

-  ! 

0  1  ...  0 

X3 

!  X2 

0  0  ...  1 

_Xn 

L^N-i 
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The  displaced  difference  of  an  arbitrary  NxN  matrix  B  is  defined  as  B-ZBZ^ .  By  means  of  the 
SVD,  the  displaced  difference  can  be  written 

S  H 

B-ZBZH  =  Xxyj 

i=l 

If  the  SVD  is  truncated  at  r  terms,  we  obtain  a  reduced  displacement  rank  r  approximation  to  B.  An 
equation  of  the  above  form  has  the  unique  solution^ 

r 

B  =1  L(Xi)U(yi) 
i=l 


where  L  is  a  lower  triangular  Toeplitz  matrix  whose  first  column  is  x^  and  U  is  an  upper  triangular 
Toeplitz  matrix  whose  first  row  is  y^ .  Then  multiplication  of  a  vector  by  B  can  be  accomplished  in 
2r  convolutions.  For  sufficiently  small  r,  a  significant  saving  of  effort  will  results  from  performing 
these  with  the  FFT,  the  number  of  operations  falling  from  OCN^),  for  the  direct  computation  of  Bg, 
to  0(2rNlog2N)  .  In  evaluating  the  product,  U(yj)g  yields  the  right  half  of  the  convolution 
between  the  top  row  of  UCyj)  and  g .  Similarly,  L(Xj)g  is  the  left  half  of  the  convolution  between 
the  bottom  row  of  Lfxj)  and  g  .  Note  that,  for  a  Toeplitz  matrix,  the  displaced  difference  contains 
non-zero  elements  only  in  the  leading  row  and  column,  is  of  rank  2  and  can  be  represented  by  the 
sum  of  two  products  of  lower  and  upper  triangular  Toeplitz  matrices. 

The  triangular  matrices  in  the  expansion  for  B  can  be  written 


Uxj)  =  [Xi|Zxi|...|ZN-ixi] 

and 

u(yi)  -  [yi|zyi|...|zN-iy|]H 

The  displaced  difference  defined  above  is  most  relevant  when  large  blocks  of  B  are  displaced  in  the 
(1,1)  direction.  In  many  cases  of  interest  (arising  in  practice  from  more  general  sampling  schemes) 
we  encounter  matrices  having  blocks  displaced  in  other  directions,  and  the  lowest  rank 
representation  will  be  obtained  when  the  displacement  reflects  the  structure  of  the  matrix.  Let  us 
define  a  generalized  displaced  difference  as 

Dp,  -  B  -  XBYH 

where  X  =  Zp  and  Y  =  Z^  . 

B  can  now  be  expanded  into  a  sum  of  products  of  lower  and  upper  triangular  matrices  (not  now  in 
general  Toeplitz): 
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and 


L(Xj)  =  [Xj  I  Xxj  I . . .  I  XN-lxj] 


U(yi)  =  [yi  I  Yyj  I . . .  I  YN-lyi]H 

DISCUSSION 

The  number  of  operations  necessary  to  form  the  product  Bg  is  approximately  2N<r+(r+l)log2N> 
where  we  have  assumed  g  is  real  and  that  the  FFT  involves  0(Nlog2N)  operations.  Note  that  for 
each  term  LUg  ,  two  FFTs  must  be  performed  since  Ug  yields  the  right  half  of  a  convolution  and 
L(Ug)  yields  the  left  half.  It  is  necessary  to  inverse  Fourier  transform  after  performing  the  first 
convolution,  select  the  appropriate  subset  of  the  data  and  then  re-transform  to  perform  the 
convolution  with  L  .  For  r  =  2  there  is  a  computational  gain  over  the  direct  matrix-vector 
multiplication  for  N  >  35  and  for  r  =  4  there  is  a  gain  for  N  >  70.  For  two-dimensional  non- 
separable  point-spread  functions,  and  r  =  4,  this  threshold  would  be  reached  by  a  9x9  image,  since 
the  imaging  operator  would  then  be  written  out  as  an  81x81  matrix.  The  computing  cost  for  a 
32x32  image  in  this  case  drops  to  about  one-tenth  that  of  the  direct  multiplication  (order  32^ 
floating-point  operations). 

Modifications  of  the  Toeplitz  imaging  matrix  A  are  induced  by  changing  the  sampling  scheme  in 
image  or  reconstruction  space;  for  example,  by  periodic  deletion  of  rows  or  columns  to  modify  the 
data  or  reconstruction  supports.  Note  that  we  are  particularly  concerned  here  with  computing  J\e 
product  Qg  ,  where  Q  is  the  regularized  pseudoinverse  of  A.  We  recall  that 
Q  =  (A^A+pi)**  AH  .  An  upper  bound  on  the  displacement  rank  of  Q  is  given  by^ 

r(Q)  <  2r(A)+3 , 

which  places  a  bound  on  the  number  of  operations  required  to  form  Qg .  We  note  also  that  for 
isoplanatic  imaging  systems  and  equal  sampling  intervals  in  data  and  reconstruction  space,  A  is 
Toeplitz  and  has  displacement  rank  2;  (A^A+pi)  can  then  be  inverted  in  0(4N2)  operations. 

It  may  not  be  convenient  from  a  hardware  point-of-view  to  take  the  data  with  the  same  spatial 
sampling  that  one  desires  to  carry  out  the  reconstruction.  If  one  forms  A  from  a  Toeplitz  matrix  T 
by  removing  alternate  rows,  then  the  sampling  interval  in  the  image  is  twice  that  in  the 
reconstruction.  For  clarity,  let  us  keep  the  same  number  of  points  in  image  and  reconstruction 
space,  so  A  is  NxN  .  Then  in  general  the  displacement  rank  of  A  can  be  as  large  as  N.  However, 
we  have  not  properly  exploited  the  structure  in  A,  which  is  along  the  direction  (1,2),  rather  than 
(1,1),  as  in  a  Toeplitz  matrix.  If  we  consider  instead  D|2(A)  =  A-ZA(Z'^)2  we  find  that  Di2(A)  is 
of  rank  <3.  In  general,  if  we  form  Aj  by  deleting  every  n***  row, 
Din(Aj)  =  A|  -2^|(Z^*)"  will  have  rank  (n+1) .  A|  can  again  be  expanded  into  a  sum  of  products 
of  lower  and  upper  triangular  matrices.  Note  that  in  this  case  the  U  matrices  are  no  longer  Toeplitz. 
However,  because  the  product  U(yi)g  yields  every  n***  point  of  the  right  half  of  the  convolution  of  yj 
and  g,  we  can  still  use  the  FFT  to  perform  the  convolution,  discarding  the  undesired  points.  If,  on 
the  other  hand,  every  n^  column  had  been  deleted  to  form  A2,  the  appropriate  displacement 


92  /  SPIC  Vot.  13S1  Digital  Image  Synthesis  and  Inverse  Optics  (1990) 


difference  would  be  D„i(A2)  =  A2  -  ,  and  a  similar  procedure  could  be  used  in  computing 

the  product  Lh  ,  where  h  =  Ug .  (The  data  points  in  h  have  also  to  be  interleaved  with  zeros  in 
appropriate  places.) 

Unfortunately,  for  A  formed  of  alternate  rows  of  a  Toeplitz  matrix  the  rank  of  D|2(Q)  and  D2i(Q) 
are  less  than  or  equal  to  the  full  rank  of  Q  even  though  Dj2(A)  =  3.  It  seems  that  one  cannot 
generally  expect  the  reduced  displacement  rank  representation  to  offer  greater  computational 
efficiency  unless  the  sampling  rate  is  the  same  in  image  and  reconstruction  spaces.  However,  in 
some  particular  cases  of  dissimilar  sampling  rates,  low  displacement  ranks  for  Q  have  been 
obtained.  For  example,  in  experiments  with  imaging  matrices  of  (sin  x/x)^  type  and  constructed 
with  the  sampling  rate  in  reconstruction  space  twice  that  in  image  space,  it  was  found  that  D2i(Q) 
was  of  much  less  than  full  rank.  The  features  of  the  matrix  which  are  responsible  for  this  property 
are  being  investigated. 


CONCLUSIONS 

A  regularized  solution  to  the  image  reconstruction  problem  has  been  discussed  and  the  use  of 
weighted  cross-validation  demonstrated  in  the  estimation  of  the  regularization  parameter  and  the 
reconstruction  support.  We  have  also  shown  how  low  displacement  rank  representations  can  be 
used  to  improve  computational  efficiency  for  imaging  matrices  of  practical  interest.  These 
techniques  are  not  limited  to  the  example  from  optical  physics  considered  here,  but  should  be 
widely  applicable  to  problems  in  the  general  area  of  information  processing. 
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IMAGE  OF  OBJECT  +  5%  GAUSSIAN  ADDITIVE  NOISE 


Figure  1:  A  scaled  image  of  an  object  17  pixels  wide.  The  image  contains  5%  Gaussian  additive  noise. 


Iog(Vs(BETASUPPORT)):  Support  Extent  Contraction;  Min(BETASUP)=(0.1,  19) 


BETA  (increasing) 

Figure  2:  Cross-validation  prediction  error  V((i,  support).  The  regularization  parameter  ranges  from  le'^to 
le*.  and  reconstruction  support  widths  range  from  25  to  1  pixel.  Tlie  minimum  occurs  for  |1  =  0.1  and 
support  width  =  19  pixels.  Tlie  object  support  was  17  pixels  wide. 
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Object  &  X-Validation  Curve  as  Support  Point  is  Excluded:  BETA=0.1 


Figure  3;  Object  and  superimposed  cross-vaIidati(xi  prediction  enor  V(P  =  0.1,  excluded_point).  Observe 
that  V  reflects  some  of  the  internal  support  stmcture  of  the  object  and  gives  a  good  estimate  of  the  object 

width. 
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Superresolution  Algorithms  for  a  Modified  Hopfield 

Neural  Network 

John  B.  Abbiss,  Bryan  J.  Brames,  and  M.  A.  Fiddy,  Member,  IEEE 


Abstract— The  purpose  of  this  paper  is  to  describe  the  imple¬ 
mentation  of  a  superresolution  (or  spectral  extrapolation)  pro¬ 
cedure  on  a  neural  network,  based  on  the  Hopfield  model.  This 
was  first  proposed  by  Abbiss  et  al.  [1].  We  show  the  computa¬ 
tional  advantages  and  disadvantages  of  such  an  approach  for 
different  coding  schemes  and  for  networks  consisting  of  very 
simple  two  state  elements  as  well  as  those  made  up  of  more 
complex  nodes  capable  of  representing  a  continuum.  With  the 
appropriate  hardware,  we  show  that  there  is  a  computational 
advantage  in  using  the  Hopfield  architecture  over  some  alter¬ 
native  methods  for  computing  the  same  solution.  We  also  dis¬ 
cuss  the  relationship  between  a  particular  mode  of  operation  of 
the  neural  network  and  the  regularized  Gerchberg-Papoulis 
algorithm. 

I.  Introduction 

There  are  several  models  of  neural  networks,  each  of 
which  has  a  structure  based  loosely  on  our  view  of 
biological  nervous  system  components  [2].  A  neural  net¬ 
work  architecture  is  one  consisting  of  a  very  large  number 
of  simple  processing  elements  densely  interconnected  by 
a  set  of  weighted  links.  Each  processing  element  updates 
its  state  by  comparing  the  sum  of  its  inputs  with  a  pre¬ 
scribed  threshold.  The  study  of  the  properties  of  neural 
networks  is  a  subject  still  somewhat  in  its  infancy,  and 
current  hardware  limitations  reduce  their  practical  im¬ 
pact.  Indeed,  it  has  been  suggested  by  Anderson  and  Ro- 
senfeld  [3]  that  they  may  not  become  useful  until  inex¬ 
pensive  special  purpose  parallel  hardware  is  available. 
Should  that  hardware  be  available,  the  question  remains 
as  to  how  one  would  make  best  use  of  a  neural  computer; 
i.e.,  how  one  should  program  or  “train”  it  to  perform  the 
tasks  required.  The  hope  is  that  some  problems  for  which 
it  is  difficult  to  find  satisfactory  algorithmic  solutions 
might  be  amenable  to  this  kind  of  computing  architecture, 
which  can  organize  itself  and  learn  what  it  is  expected  to 
accomplish. 

One  anticipated  use  of  neural  networks  is  in  autoasso- 
ciative  memory  and  in  image  (or  signal)  classification, 
recognition  or  understanding  applications;  these  are  ap¬ 
plications  that  we  believe  the  human  brain  is  particularly 
good  at  while  current  algorithms  implemented  largely  on 
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serial  machines  still  leave  much  to  be  desired.  For  most 
neural  networks,  their  learning  or  restoration  capabilities 
can  be  expressed  in  terms  of  the  minimization  of  some 
appropriate  energy  or  cost  function.  One  of  the  objectives 
of  this  paper  is  to  take  an  established  algorithm  in  image 
reconstruction  and  identify  those  aspects  of  it  that  can  be 
related  to  the  programming  requirements  that  would  be 
necessary  for  implementation  on  a  Hopfield  neural  com¬ 
puter.  From  the  analysis  of  such  a  network  applied  to  solve 
a  problem  for  which  the  cost  function  is  well  defined,  one 
might  be  able  to  assess  their  use  for  the  solution  of  a  wider 
class  of  optimization  problems. 

The  Hopfield  network  is  a  fully  connected  network  in 
the  sense  that  any  one  of  the  processing  elements  is  con¬ 
nected  to  every  other  one.  This  contrasts  with  layered  net¬ 
works,  such  as  a  multilayered  perceptron  (MLP),  in  which 
processing  elements  are  arranged  with  connections  only 
between  neighboring  layers.  This  difference  in  topology 
is  accompanied  by  differences  in  the  thresholding  func¬ 
tions  and  in  the  procedures  to  find  the  connection 
strengths.  The  Hopfield  network  operates  iteratively;  the 
connection  strengths  are  assigned  and  specify  a  cost  func¬ 
tion  which  the  iterative  procedure  minimizes.  The  MLP 
is  a  one-pass  network  once  the  connection  strengths  have 
been  “learned”  by  the  minimization  of  an  error  function 
which  quantifies  the  difference  between  the  current  and 
desired  output  states. 

It  was  pointed  out  by  Jau  et  al.  [4]  that  some  iterating 
image  restoration  processes  are  mathematically  very  sim¬ 
ilar  to  autoassociative  memory;  indeed  if  the  input  infor¬ 
mation  is  incomplete,  it  can  be  considered  as  a  key  pattern 
to  an  associative  memory.  Since  the  approach  to  image 
restoration  presented  here  was  first  proposed  [1],  there 
have  been  other  related  studies  which  we  mention  here. 
Zhou  et  al.  [S]  considered  an  energy  function  identical  to 
(6)  in  order  to  specify  network  interconnection  strengths. 
Their  application  was  the  restoration  of  grey  level  images 
degraded  by  a  shift  invariant  FIR  blur  function  and  addi¬ 
tive  noise.  Grey  level  information  was  coded  by  a  simple 
sum  of  binary  elements  and  the  network  was  serially  up¬ 
dated  with  a  stochastic  thresholding  rule  to  avoid  getting 
trapped  in  local  minima  of  the  energy  function.  Jang  et 
al.  [6]  utilized  the  optimization  properties  of  the  Hopfield 
network  in  order  to  estimate  a  matrix  inverse.  This  is 
clearly  important  in  image  restoration  problems,  as  indi¬ 
cated  by  (8).  Full  grey  level  representation  is  assumed  and 
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a  differential  mode  of  implementation  applied;  they  point 
out  that  this  method  is  similar  to  a  steepest  descent  method 
but  with  a  nonlinear  thresholding  step  at  each  iteration. 
Bai  and  Farhat  [7]  also  considered  a  cost  function  similar 
to  that  in  (6)  to  recover  images  from  limited  Fourier  data. 
In  their  approach  there  was  an  additional  constraint  on  the 
norm  of  the  derivative  of  the  estimate  as  well  as  on  the 
norm  of  the  estimate  itself.  The  increment  added  to  the 
current  estimate  was  weighted  by  a  “gain”  factor  prior 
to  thresholding,  which  was  chosen  to  ensure  that  the  net¬ 
work  energy  function  decreased  at  each  step.  Their  recon¬ 
structions  showed  advantages  for  low-signal-to-noise  ra¬ 
tio  situations  and  are  being  implemented  on  optoelec¬ 
tronic  hardware.  Winters  [8]  considers  a  norm  minimi¬ 
zation  as  expressed  by  (6)  but  without  any  explicit  regu¬ 
larization  term  included.  The  minimization  of  his  cost 
function  is  achieved  by  the  penalty  method  which  requires 
that  a  large  positive  value  is  added  to  the  cost  function, 
wherever  a  nonlinear  inequality  constraint  is  not  satisfied. 
An  adaptive  penalty  function  allows  one  to  avoid  local 
minima  and  this  complete  procedure  can  be  mapped  onto 
the  Hopfield  energy  function.  Results  showed  that  the  re¬ 
constructions  were  robust  against  noise  and  could  be  im¬ 
plemented  in  microseconds  on  an  analog  electronic  net¬ 
work,  as  compared  to  several  hours  on  a  minicomputer. 

Using  a  Hopfield  network,  our  interest  is  in  an  appli¬ 
cation  for  which  a  solution  state  evolves  through  the  min¬ 
imization  of  some  specific  cost  or  energy  function.  Once 
the  energy  function  is  defined,  one  can  determine  the  ap¬ 
propriate  connection  strengths  in  order  that  the  energy 
function  associated  with  the  network  is  the  same  as  that 
of  the  problem  under  consideration.  A  key  feature  of  a  net 
of  this  type  is  its  construction  from  a  set  of  simple  pro¬ 
cessors,  each  of  whose  states  is  determined  by  a  thresh¬ 
olding  op>eration  applied  to  a  sum  of  weighted  inputs  from 
other  processors  or  nodes.  The  properties  of  the  network 
as  a  whole  are  determined  by  the  thresholding  function 
used,  and  by  the  pattern  and  strengths  of  the  connections 
between  the  processing  elements. 


11.  The  Hopfield  Network  with  Two  State 
Elements:  Theoretical  Framework 


The  Hopfield  neural  model  [9],  [10]  allows  one  to  spec¬ 
ify  a  set  of  desired  memories  as  minima  of  a  configura¬ 
tional  energy  of  the  network.  We  assume  that  the  network 
consists  of  N  processing  elements  each  of  which  has  two 
states  and  each  of  which  has  a  thresholding  operator  that 
determines  the  states  of  the  element  from  the  total  input 
to  that  element. 

Given  an  initial  starting  configuration  or  state  of  the 
network,  each  processor  or  “neuron”  updates  its  state  ac¬ 
cording  to  a  threshold  rule  of  the  form 


N 

if  L  TiiVj  >  0  then  =  1; 

y-i 


otherwise  Vj  =  -\  { I ) 


where  the  elements  of  the  connection  matrix  T  are  formed 
according  to  the  rule 

M 

T,j  =  Z  vU'j  (ij  =  1,  2  ■  •  •  A)  (2) 

U  =  I 

and  the  cf  are  the  elements  of  the  M  memory  vectors, 
to  be  stored.  The  r’f  can  take  the  values  ±1,  and  it  is 
assumed  that  the  diagonal  term,  T,,,  is  zero  in  the  Hopfield 
model . 

This  thresholding  rule  can  be  applied  in  series  (asyn¬ 
chronously)  or  in  parallel  (synchronously).  In  the  serial 
mode,  the  rule  is  applied  sequentially  to  the  nodes  of  the 
network,  the  state  of  the  network  being  updated  after  each 
operation.  In  the  parallel  mode,  the  current  network  state 
remains  unchanged  until  the  thresholding  operation  has 
been  applied  to  every  node.  The  configurational  energy 
function  for  the  network  has  the  form  [9],  [10] 

N  N 

£  =  -  J  S  S  T.v,Vj  =  -i  v^Tv  (3) 

1  =  I  j  =  I 

where  superscript  T  denotes  transpose.  Serial  threshold¬ 
ing  will  always  minimize  this  energy  function,  provided 
that  the  7},  are  nonnegative.  If  M  is  sufficiently  small,  this 
state  will  correspond  to  the  memory  closest  in  Hamming 
distance  to  the  state  in  which  the  network  was  started. 
Parallel  thresholding  results  in  either  convergence  to  a 
stable  state  or  oscillation  between  two  states  [11]. 

This  iterative  scheme  can  be  expressed  more  concisely 
and  modified  to  permit  biasing  of  the  neuron  inputs  in  the 
following  manner.  We  let  the  state  of  the  network  after 
the  nth  iteration  be  described  by  the  N  element  vector 

where  U  is  the  threshold  operation,  T  denotes  the  matrix 
with  elements  Tij,  a  superscript  denotes  iteration  number, 
and  b  is  the  bias  vector.  The  bias  vector  incorporates 
boundary  conditions  such  as  image  data;  it  effectively 
shifts  the  decision  threshold  for  each  element.  In  this  case 
the  energy  function  minimized  by  the  network  is  of  the 
form  [10] 

E  =  -\/2  v^Tv  -  b'^v.  (4) 

The  change  in  energy  for  a  change  in  the  state  of  one 
neural  element  from  to  r*  -I-  A  I’t  is 

A£t  =  -AVkliTv  +  b)t  +  \  r«ArJ. 

Taking  to  be  zero  ensures  that  the  change  in  energy  is 
always  negative,  since  the  term  in  the  brackets  above  then 
has  the  same  sign  as  Ar’^.  If  is  nonzero,  the  term  in 
the  bracket  will  have  the  same  sign  as  A*  provided  is 
positive,  and  then  £  is  guaranteed  to  reduce;  the  conse¬ 
quences  of  varying  were  explored  in  an  earlier  publi¬ 
cation  and  verify  the  expected  behavior  [12].  If  two  (or 
more)  neurons  change  state  simultaneously,  the  change  in 
£  contains  terms  involving  products  of  the  form  -1/2 
r*,Ar/t  Af'/  (or  these  plus  higher  order  terms  if  more  neu¬ 
rons  change),  the  sign  of  which  can  vary. 
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The  two  state  representation  is  too  limited  for  a  one-to- 
one  mapping  between  elements  and  signal  or  image  sam¬ 
ples,  in  most  cases.  However,  these  simple  elements  can 
be  taken  in  groups  to  represent  grey  levels  through  a  va¬ 
riety  of  coding  schemes.  Alternatively,  either  analog  or 
more  complex  digital  processing  elements  could  be  used 
to  directly  represent  a  grey  level. 

III.  The  Superresolution  Problem 

There  are  many  applications  that  require  the  restoration 
of  a  signal  or  image  from  a  limited  discrete  data  set;  for 
example,  samples  of  the  spectrum  of  a  function,  or  of  its 
low  or  band-pass  filtered  image.  An  important  a  priori 
assumption  for  work  in  super  resolution  is  the  fact  that 
most  objects  to  be  imaged  are  of  compact  support.  This 
leads  to  the  well-known  result  that  their  spectra  are  band- 
limited  functions.  In  principle,  therefore,  one  might  hope 
to  extend  limited  spectral  data  by  means  of  analytic  con¬ 
tinuation.  This  procedure  is  notoriously  unstable  in  the 
presence  of  noise  and  does  not  provide  a  practical  solution 
to  the  problem.  One  has  infinite  freedom  in  interpolating 
and  extrapolating  limited  sampled  data;  hence,  one  is 
forced  to  approach  super  resolution  from  an  optimization 
point  of  view  [13].  The  best  that  one  can  hope  to  achieve 
is  the  specification  of  a  cost  or  energy  function  which  pos¬ 
sesses  a  unique  minimum  and  is  designed  to  incorporate 
whatever  constraints  and  a  priori  knowledge  might  be 
available  to  help  limit  the  set  of  possible  solutions  to  the 
problem,  while  retaining  desirable  and  necessaiy  solution 
characteristics.  Examples  of  constraints  include  data  con¬ 
sistency,  support  consistency  and,  perhaps,  positivity. 
The  objective  of  the  superresolution  process  is  to  obtain 
a  final  image  that  has  a  higher  spectral  or  spatial  fre¬ 
quency  content  than  the  original  data  set,  as  a  direct  con¬ 
sequence  of  incorporating  the  prior  knowledge  available 
into  the  cost  function.  It  is  a  matter  of  taste,  to  a  large 
extent,  how  one  designs  a  cost  function  in  order  to  obtain 
a  desirable  solution  to  the  problem;  i.e.,  a  superresolved 
signal  or  image  with  acceptable  properties.  The  super¬ 
resolution  problem  is  thus  transformed  into  one  of  deter¬ 
mining  the  (global)  extremum  of  a  cost  function  on  the 
assumption  that  this  solution  is  optimum. 

One  of  the  early  successes  of  a  neural  network  was  to 
find  a  good  approximation  to  the  traveling  salesman  op¬ 
timization  problem  [14].  The  superresolution  optimiza¬ 
tion  problem  can  be  mapped  onto  a  neural  network  in  two 
distinct  ways.  One  is  to  train  network  using  a  data  base 
of  superresolved  images  [15] -[17],  the  other  is  to  relate 
the  cost  function  associated  with  a  given  network  to  the 
chosen  superresolution  cost  function.  It  is  the  latter  ap¬ 
proach  that  we  adopt  here. 

IV.  Superresolution  and  Spectral  Estimation 

Most  signal  or  image  recoveiy  problems  can  be  de¬ 
scribed  by  linear  equations  of  the  form 

g(x)  =  I  A(x,  y)f(y)  dy 


where  A  is  the  system  spread  function  or  the  Fourier  trans¬ 
form  kernel,  for  example.  The  interpretation  of  the  data 
g(x)  to  obtain  information  about  the  object  /( y)  requires 
the  solution  of  a  linear  inverse  problem.  This  is  equivalent 
to  finding  the  solution  of  a  Fredholm  integral  equation  of 
the  first  kind.  It  is  well  known  that  small  fluctuations  in 
the  data  g(x)  can  lead  to  very  large  fluctuations  in  the  es¬ 
timate  of  the  unknown  function  /(>’).  This  is  a  manifes¬ 
tation  of  the  ill-posed  nature  of  the  problem  (the  inverse 
of  the  operator  A  is  not  generally  continuous)  and  some 
degree  of  regularization  is  required  in  order  to  determine 
stable  and  meaningful  solutions.  One  usually  proceeds  by 
assuming  that  the  desired  solution  belongs  to  the  space  F 
of  i^possibly  weighted)  functions,  the  regularization  re¬ 
stricting  that  solution  to  conform  to  any  a  priori  knowl¬ 
edge  available  about  the  object  whose  enhanced  image  is 
sought. 

In  practice,  the  solution  is  determined  from  a  finite  set 
of  samples  of  g(x),  and  the  data  vector  g  is  expressed  by 

g  =  Af  +  n 

where  A  is  the  imaging  operator  and  n  represents  an  ad¬ 
ditive  noise  component;  A  explicitly  contains  the  support 
constraint  of  /,  which  is  assumed  to  be  known  or  esti¬ 
mated  a  priori.  These  limited  data  can  be  regarded  as  a 
noisy  finite  set  of  bounded  linear  functionals  of  /. 

A  data-consistent  solution  exists,  however,  which  is  a 
solution  of  minimum  norm.  This  solution  is  the  data-con¬ 
sistent  g  which  minimizes  H  gP,  where  1|  1|  denotes  norm. 

The  solution  to  this  minimization  problem  can  be  writ¬ 
ten 

N 

J  =  S  (g,  t7,)M,/ai  (5) 

where  N  is  the  number  of  image  data  points  and  the  a,  ,  «, 
and  Vi  are  the  singular  values,  singular  functions,  and  sin¬ 
gular  vectors,  respectively,  pertaining  to  the  operator  A: 
Aui  =  a,  17,;  A*Vi  =  OiUi-  The  singular  values  tend  to  zero 
as  <  increases,  leading  to  an  instability  in  the  estimator.  If 
the  first  Ni  ^  N  singular  values  are  dominant,  then  the 
remainder  may  be  neglected,  but  only  at  the  expense  of 
loss  of  resolution  in 

Thus,  this  solution  is  ill-conditioned  but  stability  can 
be  restored  by  relaxing  data  consistency;  hence,  we  min¬ 
imize  the  cost  function 

£  =  11/1  r +  011  rii'-  (6) 

The  estimate  is 

N 

^  =  S  (g,  Vi)Uiai/(aj  -1-  /3)  (7) 

I  =  ! 

where  the  regularization  parameter  0  is  chosen  to  achieve 
a  compromise  between  resolution  and  stability,  and  usu¬ 
ally  requires  some  adjustment  in  order  to  establish  its  op¬ 
timal  value.  As  /3  tends  to  zero,  the  solution  becomes  more 
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data  consistent.  The  minimizer  of  this  cost  function  can 
be  computed  directly  in  matrix  form,  namely, 

J  =  [A*  A  +  &IVU*g  (8) 

where,  for  a  real-valued  matrix.  A*  becomes  A^,  and  / 
represents  the  identity  matrix. 

An  alternative  approach  to  estimating  the  object  is  to 
consider  the  minimization  of  the  cost  function  ||/  -  g'  11^ 
using  a  trigonometric  polynomial  of  the  form  [13] 

N 

g  =  2  dic<i>k 

Jfc*  I 

where  the  <t>k  form  a  basis  in  the  data  space  F  and  the 
optimal  dk  satisfy 

N 

2  [(</)„,  <t>„)F  +  =  G„  (9) 

m  *  1 

and  the  G„  are  Fourier  data  corresponding  to  the  low-pass 
filtered  image  g. 

It  is  worth  pointing  out  that  in  the  space  F  that  incor¬ 
porates  the  known  support  constraint  for  the  function  to 
be  restored,  the  three  solutions  given  by  (7)-(9)  are  equiv¬ 
alent;  expression  (9)  can  be  obtained  f^rom  expression  (8) 
[18].  Each  method  for  solution  is  more  or  less  computa¬ 
tionally  the  same  in  that  each  requires  -  0(N^)  multipli¬ 
cations;  this  was  pointed  out  earUer  in  [1]. 

We  note  that  expression  (7)  requires  on  the  order  of 
CN^  multiplications,  where  the  overhead  C  is  large  by 
comparison  with  the  other  methods.  However,  a  primary 
concern  is  the  ease  with  which  the  regularization  param¬ 
eter  0  can  be  varied;  this  can  be  done  at  the  cost  of  0(N^) 
multiplications  for  (7). 

V.  lMPLEMEt9TATION  ON  A  NEURAL  NETWORK 

We  will  now  show  how  a  superresolution  algorithm 
equivalent  to  the  previously  described  approaches  can  be 
defined  on  a  Hopfield  neural  network.  Several  issues  must 
be  addressed.  First,  it  is  necessary  to  define  the  connec¬ 
tion  matrix  from  the  cost  function.  For  some  problems 
this  cannot  be  accomplished  without  performing  more 
calculations  than  are  required  for  a  more  conventional  so¬ 
lution  to  the  problem.  The  latter  consideration  was  noted 
by  Takeda  and  Goodman  [19].  Thus,  the  computational 
load  or  complexity  must  be  considered  in  deciding  the 
merits  of  a  neural  network  solution  to  this  problem. 

Other  issues  which  must  be  addressed  center  on  modi¬ 
fications  of  Hopfield’s  formulation  to  satisfy  our  require¬ 
ments.  Hopfield  ensures  that  the  network  will  converge 
by  arbitrarily  setting  the  diagonal  of  the  connection  matrix 
to  zero.  This  is  unacceptable,  because  it  shifts  the  abso¬ 
lute  energy  minimum  of  the  network  from  the  minimum 
of  (6).  Moreover,  while  a  suitable  network  can  be  con¬ 
structed  from  two-state  elements,  one  often  requires  that 
the  reconstruction  be  represented  over  at  least  a  set  of  grey 
levels.  Thus,  it  may  be  necessary  to  combine  a  number  of 
neural  elements  to  represent  a  reconstruction  pixel.  The 
method  of  coding  the  grey  levels  depends  on  the  com¬ 


plexity  of  the  elements,  i.e.,  whether  they  can  only  take 
on  two  states,  a  bounded  continuum  of  values  such  as  0 
<  f*  <  1  (“graded  neurons”),  or  an  (effectively)  un¬ 
bounded  continuum.  Granularity  of  the  representation  will 
affect  the  convergence  properties  of  the  network;  coarsely 
quantized  systems  can  converge  to  local  energy  minima, 
yielding  less  than  optimal  reconstructions. 

We  shall  first  demonstrate  the  formal  mapping  of  a 
superresolution  algorithm  onto  a  Hopfield  network.  The 
technique  will  then  be  extended  to  fully  address  the  sec¬ 
ond  and  third  issues  on  a  two-state  network.  We  will 
briefly  examine  the  advantages  of  more  finely  quantized 
systems,  and  finally  discuss  the  relation  between  parallel 
thresholding  and  a  regularized  form  of  the  well-known 
Gerchberg-Papoulis  spectral  extrapolation  algorithm 
[20]-[22]. 

Let  us  represent  the  current  estimate  by  the  state  of  the 
network  v.  We  can  rewrite  (6)  as 

E  =  v^A^Av  -  Iv^A^g  -I-  g^g  +  0v'^v. 

Comparing  this  expression  for  E  with  that  of  the  Hopfield 
network,  (4),  gives 

T  =  -2(A^A  +  j8/)') 

b  =  2A^g  J  (10) 

where  I  is  the  identity  matrix,  and  the  g^g  term  can  be 
ignored,  since  it  represents  a  total  offset  for  E. 

Thus,  superresolution  can  be  mapped  simply  and  di¬ 
rectly  onto  a  Hopfield  network.  The  connection  matrix  is 
formed  .from  the  imaging  operator  matrix,  which  contrib¬ 
utes  information  about  the  imaging  system,  and  the 
regularization  parameter  /3,  which  sets  a  bound  on  the 
norm  of  the  final  estimate.  The  available  data  g  contribute 
only  to  the  bias  vector  b. 

For  serial  operation,  the  change  in  energy  due  to  a 
change  A  Vk  is 

A£*  =  -^Vk[{Tv  -I-  b)t  +  5  TuAt/*].  (11) 

Convergence  to  a  minimum  is  guaranteed  if  the  expres¬ 
sion  in  braces  always  has  the  same  sign  as  Ai';^.  This  can 
be  ensured  by  altering  the  diagonal  of  T to  zero;  however, 
such  a  change  is  equivalent  to  choosing  an  arbitrary  value 
for  0.  This  is  not  acceptable,. since  the  regularization  pa¬ 
rameter  should  be  chosen  to  reflect  the  noise  in  the  data 
and  to  obtain  an  optimum  reconstruction,  not  to  ensure 
convergence  of  the  algorithm. 

VI.  Superresolution  on  a  Two-State  Network 

In  this  section  we  will  modify  the  Hopfield  formulation 
so  that  the  energy  minimum  of  the  network  will  coincide 
with  that  of  (6),  while  still  decreasing  the  energy  with 
each  change  in  the  state  vector.  We  shall  find  that  this  is 
possible  by  introducing  a  two-level  threshold  in  place  of 
the  usual  single-level  one.  In  addition,  we  will  incorpo¬ 
rate  a  generalized  grey  scale  mapping  which  describes  lin¬ 
ear  transformations  of  a  state  vector  v  into  a  vector  w  hav¬ 
ing  grey  levels.  We  write  this  as  w  =  Sv,  where  S  could 
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be  a  mapping  from  an  N  element  vector,  each  of  whose 
elements  can  take  the  values  0  or  1 ,  to  an  L-element  vec¬ 
tor  (L  <  N)  whose  elements  can  take  a  wider  range  of 
values.  For  example,  if  S  represents  a  base-2  mapping, 
each  element  of  v  can  represent  a  power  of  2,  giving  a 
range  of  values  for  each  element  of  w.  The  range  of 
V  need  not  be  limited  to  {0,  1}:  we  use  this  for  the  pur¬ 
pose  of  illustration.  Other  coding  schemes  are  possible, 
such  as  clustering  or  bit-density  coding. 

The  expression  for  the  energy  is  now 

E  =  \\Aw  -  gf  +  /SlIwlP. 

Suppose  the  grey-level  vector  w  is  perturbed  by  some 
amount  Aw.  The  difference  in  energy  between  states  is 

A£  =  2Aw^{{A^A  +  0l)w  —  A^g  -t-  j  {aJa  -t-  0/)Aw} 

(12) 

or,  in  terms  of  the  neural  state  vector  v, 

A£  =  2Av^[S\A^A  +  0I)Sv  -  S^A^g 

+  ^S^(A^A  +  0I)SAv].  (13) 

It  should  be  emphasized  that  (13)  contains  no  assump¬ 
tion  about  the  range  of  values  of  v  or  of  the  updating  mode 
of  the  network.  A  restriction  to  two  states  reflects  the  de¬ 
sire  to  use  a  large  number  of  simple  binary  processing 
elements  in  neural  architectures. 

We  now  present  a  procedure  which  ensures  that  the 
change  in  energy  expressed  by  (12)  and  (13)  always  de¬ 
creases,  provid^  serial  thresholding  is  adopted.  For  a 
change  Avi^,  the  change  in  energy  AE^  of  the  network  is 
given  by  (11),  with  the  following  definitions  for  T  and  b: 

T  =  -2S\A^A  +  0I)S 

b  =  2S^A^g. 

The  grey-scale  mappings  we  are  considering  associate 
a  specific  neuron  with  one  and  only  one  image  pixel. 
Hence,  the  colunms  of  5  each  contain  only  one  element, 
and  it  is  not  difficult  to  show  that  the  diagonal  elements 
of  T  take  the  form 

Ttt  =  -1S%{A^A  -t-  mu  (14) 

where  the  nonzero  element  of  the  kth  column  of  S. 
Since  the  diagonal  elements  of  A^A  are  positive,  and  |3  is 
some  positive  quantity,  7^  is  always  negative.  Hence,  we 
can  rewrite  (13)  in  the  form 

A£*  =  -At;*{(rt;  +  b),  -  \  iralAt;*}.  (15) 

Thus,  A£t  will  be  negative  provided  {Av^f  <  A^Av^, 
where 

the  maximum  decrease  in  energy  occurring  when 

Ai;*  =  iA*.  (16) 


Thus,  we  require  that  |  Ayj|  <  [AjI  and  sgn  ( Ay^)  =  sgn 
(Ai).  For  a  binary  network,  where  t;*  e  {0,  1},  we  then 
obtain  the  following  rule  to  ensure  that  the  network  en¬ 
ergy  does  not  increase: 


r\ 

for  A*  >  1 

,Xn) 

for  |A*|  <  1 

.0 

for  A*  <  - 1 

We  consider  next  the  operation  of  a  nonbinary  network. 

VII.  Superresolution  on  a  Nonbinary  Network 

The  restriction  of  the  state  vector  v  to  binary  values 
permitted  the  simplest  possible  processing  elements  to  be 
used  in  a  neural  architecture.  With  more  complex  proces¬ 
sors  this  simple  representation  is  unnecessary  and  ineffi¬ 
cient:  the  optimal  coding  scheme  is  intimately  related  to 
the  nature  of  the  available  hardware. 

A  disadvantage  of  simple  two-level  elements  is  that  they 
can  give  a  coarsely  quantized  representation  in  recon¬ 
struction  space  which  leads  to  the  creation  of  local  energy 
minima.  There  is  still  only  one  absolute  energy  minimum, 
but  the  network  may  converge  to  a  local  minimum  of 
higher  energy.  It  should  be  recognized  that  this  behavior 
also  occurs  with  a  single  level  threshold  and  a  zero- 
diagonal  T  matrix,  and  the  resulting  reconstructions  are 
sometimes  called  “spurious  stored  states.”  A  typical  so¬ 
lution  of  this  type  is  shown  in  Fig.  1(c)  for  a  network  of 
90  two-level  elements;  a  6-b  coding  scheme  yields  15 
points  in  the  reconstruction.  Whether  such  a  reconstruc¬ 
tion  is  of  acceptable  quality  is  difficult  to  predict,  and  a 
function  of  the  needs  of  the  user.  This  difficulty  can  be 
overcome  by  using  elements  which  can  take  on  values  over 
a  continuum,  such  as  0  <  u*  <  1  (Fig.  1(e),  (f)).  These 
elements  are  similar  to  the  graded  neurons  employed  by 
Hopfield  in  a  differential  network  [10]. 

We  shall  now  examine  the  behavior  of  an  asynchronous 
network  composed  of  elements  which  can  take  on  a  con¬ 
tinuum  of  values.  Because  |  Au^il  is  no  longer  fixed  there 
is  no  need  for  a  threshold/decision  operator;  we  will  sim¬ 
ply  use  the  value  of  Avy  which  yields  the  greatest  de¬ 
crease  in  energy. 

It  was  noted  above  that  the  maximum  decrease  in  en¬ 
ergy  occurs  when 

At;*  =  A*/2. 

In  addition,  |A*|  represents  an  upper  bound  on  |  At;*|.  As 
one  approaches  the  solution,  {Tv  +  6)*  approaches  zero, 
so  this  upper  bound  decreases.  For  networks  with  a  fixed 
I  At;*|  (e.g.,  ±1  for  all  k),  one  would  expect  |A*|  even¬ 
tually  to  be  smaller  than  |  A  t;*! ,  so  no  changes  can  be  made 
to  reduce  the  energy,  even  though  the  network  is  not  yet 
at  the  global  minimum.  Thus,  the  fixed  step  methods  will 
generally  be  limited  to  some  outer  neighborhood  of  the 
global  energy  minimum  (although  one  might  arrive  at  the 
minimum). 
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Fig.  1 .  The  line  types  used  are  solid  for  an  image  or  a  neural  network  reconstruction;  dashed  for  the  object;  and  dot-dashed  for 
the  algebraic  (SVD)  reconstruction,  (a)  Object  and  incoherent  image.  Ninety  two-state  {0,  1 }  elements  are  mapped  with  a  6-b 
base  2  scheme  to  IS  pixels  in  the  object,  (b)  The  point-spread  function  of  the  imaging  system  (sinc^).  (c)  Network  using  two- 
state  neurons  converges  to  a  local  energy  minimum  after  5  cycles.  0  =  I0~’.  (d)  After  5  cycles  using  90  graded  neurons  the 
network  has  not  converged,  but  a  good  estimate  of  the  object  has  emerged,  0  =  I0~’.  (e)  Object  and  image  containing  S% 
Gaussian  additive  noise,  (f)  Neural  reconstruction  from  (e)  after  SO  cycles,  0  =  I0~*  ’.  Note  that  it  is  nearly  indistinguishable 
from  the  SVD  result. 


However,  if  we  adopt  graded  neurons  it  is  still  possible 
to  use  very  simple  elements,  yet  circumvent  the  hnite  step 
limitation.  Since  these  elements  can  take  on  a  continuum 
of  values  between  two  limits,  the  energy  is  forced  to  de¬ 
crease  at  each  step.  A  serially  threshold  netwoilc,  con- 
stracted  from  these  elements  can  therefore  reach  the  global 
energy  minimum  after  a  sufficient  number  of  cycles.  One 
could  also  dispense  with  the  coding  scheme,  and  use  a 
smaller  number  of  more  complex  elements. 

We  would  like  to  operate  the  network  in  the  synchro¬ 
nous  mode  to  make  efficient  use  of  the  network’s  paral¬ 
lelism.  If  the  Ath  neuron  changes  by 

where  convergence  is  assured  provided  that  the  conditions 
of  Section  VI  are  met. 

A  regularized  form  of  the  Gerchberg-Papoulis  algo¬ 
rithm  reads  [20]-[22] 

+  1(1  -  /3)/  - 
-  +  (7W‘">  +  b). 


Thus,  if 

Aw*  =  (Tw'"’  -I-  b)t 

parallel  operation  of  the  network  will  result  in  a  compu¬ 
tation  which  is  identical  to  the  regularized  Gerchberg-Pa¬ 
poulis  algorithm.  Since  the  latter  always  converges  [22], 
this  choice  for  the  Aw*  must  always  be  possible.  Optimal 
selection  of  the  A  w*  to  accelerate  convergence  of  the  net¬ 
work  in  the  parallel  mode  is  under  investigation. 

VIII.  Computational  Complexity  of  Neural 
Algorithm 

The  computational  complexity  associated  with  image 
reconstruction  or  superresolution  using  the  singular  value 
decomposition  of  A  to  solve  (8),  and  using  the  neural  net¬ 
work  approach,  has  been  examined  for  one-dimensional 
images  (Table  I).  The  computational  load  associated  with 
the  neural  network  is  independent  of  whether  the  thresh¬ 
olding  is  serial  or  parallel,  although  the  actual  computa¬ 
tional  time  is  obviously  less  for  parallel  thresholding.  The 
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TABLE  1 

The  Number  of  Operations  Required  for  Superresolution  by  SVD  Inversion  and  by  a  Neural 
Optimization  are  Listed  for  the  Total  Calculation,  and  for  Updating  Either  the  Image  or  the 

Regularization  Parameter 


Requirement 

Operation 

SVD 

Neural 

New  |8 

Mults 

N^  +  N 

KN^  -F  (K  -F  \)N 

Adds 

N^ 

KN^  -F  (AT  -F  l)N 

Divs 

N 

N 

New  image 

Mults 

2N^  +  N 

(K  -F  1)N^  -F  (AT  +  i)N 

Adds 

2N^ 

(AT  -F  1)N^  -F  KN 

Divs 

N 

N 

Total  operations 

Mults 

ISN^  +  3N^  +  2N 

N’  -F  (AC  -F  \)N^  -F  (AC  -F  l)N 

Adds 

N^ 

N^  +  KN^  +  KN 

Divs 

N 

N 

number  of  additions,  multiplications,  and  divisions  for 
each  technique  are  listed  in  Table  I  for  three  situations. 
The  first  row  gives  the  number  of  calculations  needed  for 
a  new  value  of  the  regularization  parameter  /3;  the  neural 
network  has  a  disadvantage  in  this  case  because  it  must 
generally  run  for  some  K  iterations.  The  second  row  gives 
the  computational  cost  involved  in  updating  the  input  im¬ 
age  data  vector  g.  The  neural  network  once  again  is  at 
somewhat  of  a  disadvantage.  However,  by  examining  the 
total  number  of  operations  from  the  beginning,  one  can 
see  that  the  neural  approach  is  substantially  more  efficient 
because  it  calculates  a  matrix  product  once  without  the 
overhead  associated  with  singular  value  decomposition; 
numerical  experiments  using  a  microcomputer  indicate 
that  the  number  of  operations  grows  as  15fV^  for  an  N 
point  image.  This  is  clearly  increasingly  significant  for 
larger  images. 

IX.  Discussion  and  Conclusions 

Neural  network  solutions  to  image  restoration  problems 
are,  therefore,  competitive  with,  but  not  necessarily  bet¬ 
ter  than,  more  traditional  methods  for  solving  the  problem 
(see  also  [22]).  We  have  shown  that  both  binary  and  non- 
binary  image  reconstruction  algorithms  can  be  imple¬ 
mented  on  very  similar  neural  architectures.  The  nonbi- 
naiy  case  can  be  based  upon  network  elements  which  take 
discrete  (typically  two-state)  or  continuous  values.  In  the 
present  application,  the  diagonal  of  the  connection  matrix 
is  nonzero  and,  for  the  discrete  element  case,  the  crucial 
thresholding  step  was  modified  in  order  to  ensure  that  the 
energy  of  the  network  decreases  at  each  step;  the  thresh¬ 
olding  step  is  unnecessary  for  the  continuous  case. 

Thus,  convergence  to  a  minimum  of  the  energy  func¬ 
tion  is  guaranteed  for  the  both  the  discrete  and  continu¬ 
ous-element  networks.  However,  this  minimum  is  unique 
only  in  the  continuous  case,  since  discretization  of  the 
element  states  introduces  local  minima.  Convergence  has 
been  demonstrated  in  this  paper  only  for  the  case  of  serial 
thresholding;  for  parallel  thresholding  in  the  continuous 
case,  the  computation  reduces  to  the  regularized  Gerch- 
beig-PapouIis  algorithm  for  a  particular  choice  of  the  in¬ 
crements  in  the  states  of  the  network  elements. 
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ABSTRACT 

We  present  an  analysis  and  computational  results  relating  to  the  regularized  restoration  of  subpixel  information  from 
undersampled  data.  Tlie  method  makes  use  of  a  small  set  of  images  in  various  stages  of  defocus.  An  iterative  implementation 
permits  the  incorporation  of  a  non-negativity  constraint  The  problem  we  consider  is  fundamentally  under-determined,  but 
useful  results  can  be  obtained  in  reasonably  low  noise  conditions. 

1.  INTRODUCTION 

The  investigations  discussed  here  form  part  of  a  program  whose  subject  is  the  enhancement  of  images  obtained  from  space- 
based  remote  sensors.  For  the  present  purpose,  these  images  are  assumed  to  consist  of  quantized  data  from  a  fixed  two- 
dimensional  set  of  sensors,  such  as  a  CCD  array.  Typically  for  these  arrays,  the  Airy  disc  is  smaller  than  one  pixel.  Thus,  for 
reasonably  fast  and  well-corrected  optics,  the  conventional  limit  on  system  resolution  is  likely  to  be  the  result  of  the  spatial 
sampling  associated  with  the  pixel  size;  i.e.,  the  integration  of  the  light  energy  in  the  image  over  the  area  represented  by  each 
pixel. 

For  clarity,  we  consider  the  case  of  one-dimensional  imaging  with  an  incoherent  source.  Let  the  system  point  spread  function 
(psf)  be  represented  by  the  continuous  imaging  operator  L  .  Then,  for  an  isoplaoalic  system  at  Uie  absence  of  noise,  the 
image  g  of  an  object  f  is  given  by  the  convolution 

g(y)  =  Js  L(x  -  y)  f(x)  dx  (1) 

where  S  is  the  support  of  f  .  The  output  of  the  k*  detector  (pixel),  extending  from  to  Y|j+j,is 

fYk+l 

gk  =  JYk  8(y)dy 


-  L(x  -  y)  f(x)  dx 

Interchanging  the  order  of  integration  (which  is  certainly  permissible  with  the  physical  functions  considered  here)  and  making 
the  definition 

rYk+l 

ak(x)  =  Jy,^  L(x-y)dy  , 

we  obtain  for  the  pixel-integrated  image 

gk  "•  /s  M*)  fW  dx,  k  “  1, 2,  •  K  (2) 
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We  now  discretize  the  object  into  a  set  {fj,  j  =  1.2  J}  over  equal  intervals  and  write,  as  an  approximation, 

j 

gk  =  atjfj,  k=  1.2,  -.K 

where  a|(j  is  the  integral  of  a^Cx)  over  the  j***  interval.  In  matrix  form,  equation  (3)  can  be  written 

g  =  Af 


(3) 


where  A  is  a  KxJ  matrix. 

Since  the  a|g-  can  be  determined  from  the  system  psf  and  the  pixel  array  structure,  the  deconvolution  problem  represented  by 
equation  (3)  can,  in  principle,  be  solved  and  the  set  {fj}  leconstracted  approximately  from  an  equal  set  of  measurements  of 
the  gk. 


In  practice,  the  restoration  problem  is  made  much  more  difficult  because  of  background  and  intrinsic  noise,  array 
imperfections  and  errors  or  uncertainties  in  the  knowledge  of  the  optical  properties  of  the  system.  When  pixel  integration  is 
over  a  significant  part  of  the  system  psf,  achieving  even  a  modest  degree  of  enhancement,  using  a  single  image,  poses 
intractable  difnculties  in  the  presence  of  these  perturbations.  However,  the  necessary  additional  information  can  be  derived 
from  multiple  differing  images  of  the  given  object  After  a  brief  discussion  of  the  characteristics  of  the  inverse  problem 
represented  by  equation  (3),  we  shall  give  a  detailed  description  of  a  specific  method  for  acquiring  this  information. 

2.  REGULARIZED  IMAGE  RESTORATION 

There  are  many  applications  that  require  the  restoration  of  a  signal  or  image  from  a  limited  discrete  data  set;  for  example, 
samples  of  the  spatial  or  temporal  spectrum,  or  of  the  object's  low  or  band-pass  tillered  image.  An  important  a  priori 
assuiiq)tion  in  image  restwation  is  that  the  object  or  objects  are  of  compact  support  This  leads  to  the  well-known  result  that 
their  spectra  are  bandlimited  functions.  In  principle,  therefore,  one  might  hope  to  extend  limited  spectral  data  by  means  of 
analytic  continuation,  and  then,  by  Fourier  transfoimation,  obtain  an  enhanced  image.  This  procedure  is  notoriously  unstable 
in  the  presence  of  noise  and  does  not  provide  a  practical  solution  to  the  problem. 

This  central  difticulty  can  be  expressed  in  another  way.  Most  signal  or  image  recovery  problems  can  be  described  by  linear 
equations  of  the  form 


g(x)  =  /  A(x,y)  f(y)  dy 

where  A  is  the  system  point  spread  function  or  the  Fourier  transform  kernel,  for  example.  The  inteipretation  of  the  data  g(x) 
to  obtain  information  about  the  object  f(y)  requires  the  solution  of  a  linear  inverse  problem.  This  is  equivalent  to  finding  the 
solution  of  a  Fredholm  integral  equation  of  the  first  kind.  It  is  well-known  in  such  a  case  that  small  fluctuations  in  the  data 
g(x)  can  lead  to  very  large  fluctuations  in  the  estimate  of  the  unknown  function  f(y).  This  is  a  manifestation  of  the  ill-posed 
nature  of  the  problem  (the  inverse  operator  is  unbounded)  and  some  method  of  stabilization  is  needed  to  determine  useful 
solutions. 

The  problem  of  a  lack  of  continuous  dependence  on  the  data  can  be  overcome  by  one  of  the  various  techniques  of 
regularization  [1,  2].  Essentially,  the  ill-posed  problem  is  replaced  by  a  related  well-posed  one.  chosen  to  be  i^ysically 
meaningful  and  to  possess  the  necessary  prc^rties  of  convergence  and  stability.  Thus  we  may  change  the  concept  of  a 
solution,  or  the  Hilbert  spaces  of  which  tire  object  and  image  are  elements,  m  their  topologies,  or  the  operator  itself.  The 
tedmique  we  shall  use  belongs  to  the  last  category. 

We  impose  physically  reasonable  constraints  on  the  permitted  solutions.  If  e  is  a  measure  in  the  norm  sense  of  the  noise  in 
the  unage  (noim  in  tte  appropriate  Hilbert  space  is  denoted  by  II  •  1 1  h  )  and  if  C  is  some  constraint  operator  with  E  a 
known  bound,  we  shall  require  that  all  possible  reconstructions  f  satisfy 
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||g-A/’llG<e  and  11  C/-  11f<E 

C  may  be  used,  for  example,  to  impose  smoothness  on  the  reconstruction  or  to  weight  the  reconstruction  support.  If  C  is  the 
identity  operator,  E  is  a  bound  on  the  norm  of  the  reconstruction.  We  combine  the  constraints  quadratically  and  minimize  the 
functional 


llg-A/’llc  +  pll  Cfllp 


where  P  =  e^/E^.  Note  that  smaller  values  of  P  are  equivalent  to  demanding  greater  fidelity  between  the  reconstruction  and  the 
data;  greater  values  place  more  emphasis  on  the  property  of  the  reconstruction  associated  with  C.  In  the  present  discussion  we 
shall  take  C  =  I.  Theminimizer  fp  can  then  be  expressed  in  either  of  the  forms 


/p  =  (A*A  +  pI)->  A*g 
or 

/p  =  A*  (AA*  +  pI)-«  g 


(4) 


where  A*  is  the  operator  adjoint  to  A.  The  inverses  of  the  bracketed  operators  will  always  exist,  since  the  eigenvalues  of  the 
symmetric  operators  A* A  and  AA*  are  non-negative.  The  operator,  the  image  and  the  reconstruction  will  consist  in  practice 
of  finite  arrays  and  equation  (4)  becomes  a  matrix  equation. 


3.  SUB-PIXEL  RESOLUTION  FROM  MULTIH^E  DEFOCUSSED  IMAGES 


We  can  state  the  deconvolution  problem  in  the  more  general  case  where  we  are  given  a  set  of  R  differing  noisy  images  of  the 
same  object  over  the  same  pixel  array. 

Then  the  image  consists  of  the  vector  of  pixel  outputs  given  by  the  equation 

g«  =  A,f+nfr)  (5) 

where  nl*^)  represents  some  additive  noise  vector.  (Other  forms  of  noise  can  be  accommodated  by  appropriate 
modifications  of  equations  (4)  and  (S),  but  we  shall  not  address  this  question  here.)  We  are  required  to  estimate  f  from  the 
set  {gfr>,  r-1,2, 

We  obtain  the  solution  by  assembling  the  image  set  into  one  ctnnposite  image,  and  the  corresponding  matrices  of  integrated  psf 
samples  a<j|[)  into  a  single  imaging  matrix.  The  regularized  solution  is  then  again  given  by  equation  (4).  The  appropriate  value 
for  the  regularization  parameter  p,  which  is  closely  related  to  the  signal-to-noise  ratio  in  the  data,  can  be  estimate  in  several 
ways;  for  exanqile,  by  the  method  of  weighted  cross-validation  [3]. 

The  image  set  required  to  achieve  sub-pixel  resolution  can  be  derived  by  various  means.  Stark  and  Oskoui  [4]  discuss  an 
object  reconstruction  technique  which  uses  a  set  of  images  differing  from  one  another  by  rotation  or  lateral  translation  of  the 
pixel  array.  It  is  evident  that  acquiring  a  sequence  of  data  sets  by  lateral  displacement  of  the  detector  (or  image)  through  some 
fraction  of  a  pixel  at  each  step  allows  one  to  sample  the  image  as  finely  as  is  desired.  (For  the  two-dimensional  image,  the 
di^lacements  can  be  in  any  direction.)  Rotational  displacement  similarly  permits  arbitrarily  fm:  sampling.  They  consider  the 
specific  case  in  which  the  system  psf  is  much  smaller  than  a  pixel,  so  that  resolution  in  the  data  is  governed  almost  conqrletely 
by  pixel  integration  rather  t^  the  optical  properties  of  the  system. 

An  alternative  method,  which  we  shall  adopt,  is  to  use  an  image  set  obtained  with  differing  point-spread  functions.  These 
could  be  generated  by  separate  optical  systems,  or  conceivably  at  different  wavelengths.  The  system  psf  can  also  be 
convenienUy  altered  by  varying  the  degree  of  defocus;  implementations  of  this  method  could  depend  on  a  single  detector 
translated  into  the  chosen  planes  of  defocus,  or  a  system  of  beamsplitters  and  detector  arrays  in  appropriate  locations. 

If  the  reconstruction  procedart  is  to  be  effective,  the  various  images  must  contain  significantly  different  information,  which 
implies  that  the  point  spread  functions  must  differ  appreciably  over  the  scale  of  a  pixel.  We  shall  assume  that  the  images  are 
formed  on  the  same  array,  or  arrays  with  identical  characteristics,  and  that  the  object  field  is  spatially  and,  if  necessary. 
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temporally  invariant  from  one  image  to  another.  We  shall  also  assume  that  the  point  spread  functions  are  accurately  known 
and,  for  the  purposes  of  the  algorithm  used  later,  that  the  images  are  formed  from  incoherent  radiation.  We  do  not  require  that 
any  of  the  point  spread  functions  are  much  smaller  than  a  pixel;  for  the  illustration  presented  below,  the  psf  at  focus  is  about 
half  the  size  of  a  pixel.  In  a  practical  implementation,  the  images  would  not  be  centered  on  the  pixel  arrays  in  exactly  the  same 
way;  i.e.,  there  would  be  some  lateral  translation  between  the  images.  The  effect  of  including  lateral  translations  has  not  yet 
been  investigated,  but  one  would  not  expect  reconstruction  quality  to  suffer  under  these  circumstances. 

4.  ILLUSTRATION 

To  illustrate  these  ideas,  we  include  some  results  from  a  numerical  simulation.  The  object  field  consists  of  a  group  of  delta- 
function-like  incoherent  radiators.  We  shall  show  that  from  a  small  set  of  images,  one  at  focus,  the  others  at  various  stages  of 
defocus,  object  locations  can  be  well  recovered,  even  when  several  of  the  objects'  geometrical  images  lie  within  a  single  pixel. 
The  algorithm  is  based  on  the  regularized  formula  of  equation  (4);  in  addition,  the  calculation  is  iterated  a  small  number  of 
times  to  enforce  non-negativity  on  the  reconstruction. 

The  reconstruction  is  made  initially  into  a  spatial  region  defined  by  the  central  lobe  of  the  focussed  image,  but  using  a  finer 
grid.  No  prior  assumptions  are  made  about  the  locations  or  the  number  of  objects  within  this  region.  For  small,  relatively 
isolated  sources,  some  ringing  will  occur  in  the  reconstmction,  with  associated  negative  pixel  values.  The  support  is 
progressively  reflned  by  eliminating  these  pixels  at  each  iteration  until  an  entirely  positive  reconstruction  is  obtained.  The 
smallest  object  space  which  is  consistent  with  the  image  data  will  yield  the  best  reconstmction,  and  the  problem,  initially 
underdetermined,  becomes  finally  an  overdetermined  one. 

A  modified  form  of  this  scheme,  which  would  be  appropriate  for  more  extended  objects,  includes  a  weight  matrix  which 
biasses  the  next  iteration  against  those  pixels  with  negative  values.  The  computation  is  significantly  slower  in  this  case,  since 
the  size  of  the  reconstruction  space  remains  constant.  We  also  note  the  possibility  of  using  a  regularized  form  of  the  non¬ 
negative  least-squares  algorithm  of  Lawson  and  Hanson  [5].  The  relative  performance  of  this  procedure,  also  of  course 
iterative,  has  not  been  fully  evaluated. 

In  this  example,  the  object  field  consists  of  eleven  highly  localized  sources  of  equal  intensity,  distributed  over  a  3x3  block  of 
image  pixels;  see  Figure  1.  (It  should  be  noted  that  the  reconstmction  grid  used  did  not  coincide  with  the  object  grid;  hence  the 
objects  cannot  appear  as  single-pixel  "points"  in  the  reconstmction.)  The  central  lobe  of  the  system  psf  was  about  half  the 
width  of  a  pixel.  Four  independent  images  containing  equal  energy  were  generated  over  a  7x7  block  of  pixels,  the  first 
corresponding  to  a  focussed  system,  the  others  at  various  stages  of  defocus.  A  method  was  devised  for  choosing  defocus 
conditions  with  significantly  different  information  content  Starting  with  the  focussed  image,  with  associated  imaging  matrix 
Ag,  we  wish  to  find  a  defocussed  image  whose  imaging  matrix  Aj  is  as  independent  as  possible  from  Ag.  The  criterion  used 
for  this  purpose  was  the  magnitude  of  the  condition  number  (the  ratio  of  the  largest  to  the  smallest  singular  values)  of  the 
matrix  formed  from  the  center  columns  of  Ag  and  Aj.  The  range  of  defocus  over  which  the  search  was  made  was  from  0  to  4 
waves.  The  combination  selected  was  that  possessing  the  smallest  condition  number.  Knowing  Ag  and  A),  A2  was 
determined,  and  finally  A3.  The  merit  functions  (reciprocals  of  condition  numbers)  calculated  for  combinations  of  two,  three 
and  four  defocus  levels  are  shown  in  Figure  2.  Note  that  a  zero  occurs  whenever  the  variable  degree  of  defocus  coincides  with 
that  corresponding  to  one  of  the  other  imaging  matrices;  then  the  composite  matrix  is  singular  and  its  condition  number 
bcomes  unbounded.  The  degrees  of  defocus  chosen  for  the  four  images  in  this  case  were  of  magnitude  0,  0.8,  1  6  and  2.4 
waves.  The  corresponding  images  are  shown  in  Figures  3-6.  The  only  readily-identifiable  feature  in  the  focussed  image  is  that 
the  central  pixel  is  brighter  than  the  surrounding  ones. 

Reconstructions  were  performed  over  the  region  defined  by  the  central  block  of  3x3  pixels,  with  sampling  seven  times  as  fine 
as  that  in  the  data.  Thus  the  reconstruction  space  initially  consisted  of  441  points,  while  the  four  images  provided  a  total  of 
1%  data  points.  The  reconstmctions  obtained  when  each  of  the  images  was  corrupted  by  additive  Gaussian  noise  with 
standard  cteviation  equal  to  1%  of  the  mean  pixel  content  is  shown  in  Figure  7.  All  of  the  objects  are  located  close  to  their  true 
subpixel  positions.  Figure  8  shows  the  result  obtained  when  the  noise  level  is  increased  to  5%.  The  reconstmction  is  still 
generally  accurate,  although  there  is  now  some  distortion  in  object  location  and  one  or  two  small  artifacts  have  begun  to 
appear.  A  signal-to-noise  ratio  can  be  defined  as  the  ratio  for  e^cb  image  of  the  sum,  on  a  pixel-by-pixel  basis,  of  the  signal 
power  to  the  sum  of  the  noise  power.  At  the  5%  level,  this  quantity  varied  between  33  and  28  dB.  This  example  was  designed 
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to  be  reasonably  challenging,  and  spreading  the  objects  further  apart,  or  reducing  their  number,  considerably  increases  the 
algorithm's  robustness  against  noise. 


5.  CONCLUSIONS 

A  method  for  recovering  detail  at  the  sub-pixel  level  from  a  small  set  of  images  in  various  stages  of  defocus  has  been  discussed 
and  demonstrated,  using  an  algorithm  which  incorporates  regularization  to  counter  the  destabilizing  effects  of  noise.  The 
iterative  version  of  this  algorithm,  which  permits  the  inclusion  of  a  non-negativity  constraint,  is  particularly  effective  at 
recovering  accurate  object  support  estimates.  It  is  therefore  an  appropriate  technique  to  use  when  it  is  known,  a  priori,  that  the 
object  field  consists  largely  of  small  well-separated  targets.  The  sensitivity  to  noise  of  the  method  deserves  more  detailed 
investigation.  A  quantitative  comparison  of  its  performance  with  methods  which  make  use  of  laterally  translated  or  rotated 
image  sets  would  also  be  of  considerable  interest  In  practice  there  would  inevitably  be  some  lateral  shift  between  images, 
even  if  they  are  formed  on  the  same  array,  and  a  hybrid  scheme  might  prove  to  be  the  most  robust  in  the  presence  of  noise. 
Ultimately,  of  course,  the  performance  of  any  restoration  algorithm  must  be  limited  by  the  information  content  of  the  image 
set  What  should  be  sought,  therefore,  is  the  eiKoding  scheme  which  most  efficiently  exploits  the  total  information  carried  by 
the  incident  radiation. 
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Fig.  1. 
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Fig.  2.  Merit  functions 
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Fig.  3.  Image  at  0.0  waves  defocus 


Fig.  4.  Image  at  0.8  waves  defocus 
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Fig.  5.  Image  at  1.6  waves  defocus 


Fig.  6.  Image  at  2.4  waves  defocus 
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Fig.  8.  Reconstruction  from  image  with  5% 
added  noise  (beta  =  10'^) 
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ABSTRACT 

A  new  algorithm  for  the  restoration  of  extended  images,  the  Regularized  Pseudoinverse  Deconvolution  (RAPID) 
algorithm,  is  proposed.  The  algorithm  consists  of  expanding  the  regularized  pseudoinverse  of  the  imaging  operator 
into  a  sequence  of  terms  which  can  be  easily  implemented  using  Fourier  processing  techniques.  The  first  term  of  the 
expansion  is  closely  related  to  generalized  V^ener  filtering  if  the  point  spread  function  is  shift-invariant.  The  other 
terms  in  the  expansion  are  correction  terms  which  cire  small  when  the  point  spread  function  is  shift-invariant,  as  is  the 
case  with  many  imaging  systems.  Even  when  the  point  spread  function  of  the  imaging  system  is  space-variant,  such  as 
with  a  partially  obscured  imaging  system  or  a  system  with  severe  aberrations,  the  correction  terms  are  both  few  in 
number  and  easily  implemented. 


1.  INTRODUCTION 

In  the  absence  of  any  other  degrading  effects,  the  performance  of  an  optical  system  is  ultimately  restricted  by 
diffraction.  The  finite  extent  of  the  entrance  pupil  imposes  a  fundamental  upper  limit  on  the  system's  spatial  frequency 
response.  The  image  quality  of  most  operational  s)^tems  will  not,  however,  approach  this  theoretical  limit  very  closely. 
It  is  possible  that  the  design  or  construction  will  be  flawed,  as  in  the  case  of  the  Hubble  Space  Telescope,  through 
defective  manufacture,  assembly  or  quality  assurance  procedures.  In  addition,  aging  of  components  will  almost 
certainly  compromise  sensor  performance  at  some  stage  in  its  lifetime,  and,  as  in  spacebome  operation,  replacement 
may  not  be  a  simple  task.  The  detector  itself  may  impose  limitations;  for  example,  where  a  CCD  array  is  used, 
information  is  lost  in  the  inter-pixel  areas,  and  image  energy  is  integrated  over  the  active  area  of  each  pixel.  Other 
degrading  factors  will  include  defective  pixels,  noise  in  the  CCD  array  and  electronic  subsystems,  and  a  possibly 
obtrusive  background.  The  methods  of  image  restoration  considered  here  were  originally  aimed  at  achieving 
performance  beyond  the  diffraction  limit  ^  but  are  in  fact  capable  of  compensating  simultaneously  or  separately  for 
aberrations  induced  by  the  optical  components  and  for  the  limitations  of  the  detector.  They  arc  inherently  robust  and 
possess  valuable  noise-suppressing  properties. 
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The  assumption  is  made  that  the  overall  effect  of  the  optics  can  be  described  as  a  possibly  time-  and  shift-variant 
blurring  of  the  image  due  to  diffraction  and  aberrations.  Thus,  at  any  instant  of  time,  the  point  spread  function  may 
change  across  the  sensor  field-of-view.  It  can  be  assumed,  however,  that  at  any  given  point  in  the  image  the  point 
spread  function  is  effectively  determined  by  its  spatial  location  and  is  time-invariant  over  the  integration  time  of  the 
sensor  and  then  undergoes  a  change  at  a  later  time.  This  is  the  case  with  the  Hubble  Space  Telescope.  It  will  be 
assumed  that  the  set  of  point  spread  functions  is  known  or  can  be  measured.  For  the  Hubble  wide-field  planetary 
camera  there  are  two  primary  components  to  the  point  spread  function;  one  due  to  diffraction  by  the  aperture 
obstructions  and  the  other  due  to  spherical  aberration.  The  point  spread  function  is  observed  to  be  locally  shift- 
invariant,  and  the  image  can  be  considered  to  be  created  by  the  summation  across  the  entire  field-of-view  of  the 
segmented  set  of  localized  point  spread  functions  convolved  with  objects  in  the  corresponding  parts  of  the  field. 

However,  for  a  simple  refracting  telescope  designed  for  observation  of  extended  objects,  the  aberrations  causing 
image  degradation  can  be  highly  space-variant.  For  an  extended  object,  image  segmentation  introduces  undesirable 
edge  effects  at  the  block  boundaries  when  locally  shift-invariant  approximations  of  the  globally  shift-variant  point 
spread  function  are  used  to  process  the  separate  blocks.  In  addition,  post  processing  interpolation  or  iteration  is 
required  to  smooth  the  block  boundaries  in  the  final  composite  image  In  the  approach  proposed  in  this  paper,  the 
pxjint  spread  function  is  allowed  to  vary  continuously  over  the  whole  image  and  a  special  decomposition  of  the  image 
reconstruction  operator  is  performed  which  permits  Fourier  techniques  to  be  used  to  process  the  entire  image. 

The  overall  optical  system  imaging  equation  can  be  written  as  a  Fredholm  integral  equation  of  the  first  kind. 
Solutions  to  ill-posed  problems  of  this  type  are  known  to  be  numerically  unstable  Additionally,  it  is  anticipated  that 
the  image  will  be  spatially  sampled  by  a  solid-state  sensor  which  will  introduce  spatial  integration,  discretization  and 
associated  noise  processes.  Thus,  after  scanning  the  image  into  a  vector,  the  integral  representing  the  continuous  image 
can  be  rewritten  as  a  matrix  expression.  In  general,  the  presence  of  the  sensor  noise  takes  the  measured  image  vector 
out  of  the  span  of  the  columns  of  the  kernel  matrix,  which  is  typically  highly  ill-conditioned.  Thus,  when  it  is  desired  to 
estimate  what  the  image  would  have  been  in  the  absence  of  aberrations  or  with  less  diffraction  than  the  instrument  can 
provide,  techniques  derived  from  regularization  theory  are  required  to  restore  stability  to  the  reconstruction.  By 
introducing  a  suitable  error  criterion  (based  on,  e.g.,  a  vector  norm),  images  can  be  constructed  which  are,  in  terms  of 
the  chosen  criterion,  closer  to  the  undistorted  geometrical  image  of  the  object  than  the  detected  image  data. 

To  the  extent  permitted  by  the  noise  in  the  image,  in-band  effects  can  usually  be  removed  by  some  form  of 
pseudoinverse  filter  However,  detector  pixellation  and  the  finite  aperture  of  any  system  set  resolution  limits  not  so 
easily  overcome,  and  a  method  for  achieving  spectral  extrapolation  has  to  be  devised.  The  spatial  spectrum  of  the 
object  is  the  Fourier  transform  of  its  amplitude,  in  the  coherent  case,  or  its  intensity,  in  the  incoherent  case.  If  the  object 
is  known  to  be  of  finite  extent,  its  Fourier  transform  is  an  analytic  function,  and  the  out-of-band  part  of  the  spectrum 
can  in  principle  be  fully  recovered  by  analytic  continuation  ^  of  the  image  spjectrum  after  removal  of  any  in-band 
distortion.  The  inverse  Fourier  transform  of  this  extended  spectrum  would  then  yield  a  perfect  image  of  the  original 
object.  Equivalently,  one  could  attempt  to  solve  directly  the  equation  describing  the  imaging  process.  This,  however, 
involves  the  inversion  of  an  ill-conditioned  matrix,  and  the  restoration  process  is  intrinsically  unstable,  even  small 
amounts  of  noise  rendering  the  results  meaningless.  These  difficulties  may  be  surmounted  by  applying  the  methods  of 
regularization  theory  developed  to  deal  with  ill-posed  problems  of  this  type;  the  solution  is  derived  by  means  of  a 
constrained  least-squares  procedure  in  which  a  regularization  parameter  plays  an  essential  role.  Stability  in  the 
restored  image,  which  is  computed  via  the  regularized  pseudoinverse  of  the  imaging  matrix,  is  controlled  by  this 
parameter.  Its  optinuil  value  depends  on  the  signal-to-noise  ratio  in  the  data. 
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2.  NATURE  OF  THE  PROBLEM 


We  wish  to  estimate  the  object  f  from  an  image  g,  given  that 


g  =  Af+r 


where  A  is  the  imaging  operator  and  r  represents  the  corrupting  effect  of  additive  noise.  For  clarity  in  the  analysis, 
we  consider  the  one-dimensional  case  with  the  operator  A  given  in  integral  form  by 

(A/)  (y)  =  J  A(x,y)f(x)dx,  c<y<d.  (2) 

In  the  absence  of  noise,  Eq.(l)  becomes  a  Fredholm  equation  of  the  first  kind,  in  which  the  unknown  function  appears 
only  under  the  integral  sign.  We  can  identify  the  sources  of  difficulty  in  solving  this  equation  in  the  presence  of  noise  or 
other  jjerturbations  (such  as  computer  round-off  error)  by  means  of  a  singular  function  analysis  ®. 

We  expand  the  kernel  of  the  integral  in  terms  of  the  singular  functions  u.(x)  and  (y) ,  orthonormal  systems  in 
object  and  image  spaces  respectively,  and  the  singular  values  a.: 


A  (X,  y)  =  ^  a.  (X)  V.  (y) . 


The  object  and  image  can  be  expanded  in  the  forms 


fix)  =  2^/,  K..I 


s(y)  =  !>,•(») 


where  the  coefficients  are  related  to /(x)  and  g(y)  by  the  integral  formulae 


/,•  =  J*  f{x)u.{x)dx 

=  I  giy)v^iy)dy. 


In  the  noiseless  case  (r  =  0) ,  we  find  from  Eqs.  (2)  and  (6) 


(Af)  (y)  =  /.  v.{y). 


3 


It  follows,  using  Eqs.(5)  and  (8),  that 
and  hence 


«■  =  o.  {■ 


(9) 


(10) 


i  =  l  * 

Thus  the  object  function  can  in  principle  be  perfectly  reconstructed  from  the  set  }  of  image  coefficients. 

Now  consider  the  effects  of  noise.  By  expanding  r  (y)  in  terms  of  the  (y)  ,  we  can  derive  the  contribution  of 
the  noise  to  the  new  image  coefficients: 


The  estimate  of  }{x)  is  now 


r'.  =  a.  f.  +  r.. 
*  z  t  n  i 


r. 


fix)  =  fix)  +  V  ^  u.ix). 


Z  =  1 


(11) 


(12) 


Image  formation  is  characteristically  described  by  an  integral  transform  of  convolution  type,  i.e.,  is  a 
convolution  operator.  Its  singular-value  spectrum  typically  decays  asymptotically  at  an  exponential  rate  Since  the  r. 
will  in  general  decrease  less  quickly,  the  sum  in  Eq.(12)  will  be  divergent  and  no  bound  will  exist  for  the  'distance'  (in 
the  sense  of  some  appropriate  metric)  between  the  true  object  and  the  reconstruction.  The  effect  of  the  noise  on  the 
reconstructed  image  is  a  manifestation  of  the  fact  that  convolution  is  a  strongly  smoothing  process  -  closely  similar 
images  can  correspond  to  widely  differing  objects.  Thus  image  restoration  is  an  ill-posed  problem,  small  perturbations 
in  the  data  causing  laige  changes  in  the  solution  represented  by  Eq.(12). 


3.  THE  REGULARIZED  SOLUTION 


The  methods  of  regularization  theoiy  ^®can  be  exploited  to  convert  the  problem  to  a  related  well-posed  one,  i.e., 
one  for  which  the  solution  exists,  is  unique  and  depends  continuously  on  the  data.  Since  we  shall  later  be  concerned 
with  computed  reconstructions,  we  henceforth  consider  the  problem  in  its  finite  discretized  form;  thus  the  imaging 
operator  A  becomes  a  matrix.  Although  strictly  we  should  introduce  new  symbols,  for  convenience  we  continue  to 
use  /,  g,  and  A  to  denote  the  discrete  forms  of  object,  image,  and  imaging  operator.  Generally  feC^  ,  g  €  C”*  and 


To  regularize  the  problem,  we  shall  modify  A.  We  impose  constraints  on  possible  solutions  /  'by  requiring  that 

\\Af'-g\\^<t^  (13) 


where  e  is  some  suitable  measure  of  the  noise  in  the  image,  and  that 


ll/'II^S^^  (14) 
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where  ^  is  some  suitable  measure  of  the  permitted  'signal  strength'  of  the  solution.  (  |1  •  ||  denotes  norm  in  the 
Hilbert  spaces  associated  with  object  and  image.)  We  combine  these  constraints  and  minimize 


\\Af-g\\^+M  r  11^ 


where  the  regularization  parameter  P  is  given  by 


The  minimum-norm  solution  to  this  constrained  least-squares  problem  is  given  by 


where 


/p  -  ^p« 


+  A^. 


We  note  the  relatioirship  of  (which  we  shall  call  the  regularized  pseudoinverse)  to  A  ,  the  Moore-Penrose 
pseudoinverse 

A'*’  =  lim  At. 


tirn  « 

p^o  P 


H  H 

The  inverse  of  {A  A  +  pi)  always  exists,  since  A  A  is  non-negative  definite  and  the  regularization  parameter  p 
is  positive.  A  value  of  P  should  be  chosen  which  balances  data  fidelity  against  smoothness  in  the  reconstruction. 
Methods  are  also  available  for  determining  P  from  the  image  data  themselves 

4.  SINGULAR  VALUE  AND  FOURIER  DECOMPOSITIONS 
It  is  often  convenient  to  compute  the  regularized  pseudoinverse  via  the  singular  value  decomposition  (SVD)  of  A: 


where 


A  =  UZV^ 


U^U  =  I  ,  V^V  =  =  I 

m'  n 


Z  =  diag  (Cj,  <J^, o^) ,  a.  >  0. 


The  singular  values  are  assumed  to  have  been  arranged  in  descending  order  of  magnitude 

G-  >a„  >a,  > ...  >a  . 

12  3  n 


Then  we  find 


/p  =  V  Eg  U”  « 
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where 


=  diag 


a^  +  P 

I  ^ 


This  representation  is  useful  when  the  behavior  of  the  reconstruction  as  a  function  of  the  regularization  parameter  is 
being  studied.  Regularization  can  equivalently  be  achieved  by  setting  p  to  zero,  and  simply  truncating  the  singular 
value  series  at  a  point  which  is  dependent  on  the  noise  level 

An  ah  initio  computation  of  /p  via  Eq.(23)  requires  the  SVD  of  A  followed  by  two  matrix-vector  multiplications.  If 
the  regularized  pseudoinverse  can  be  precomputed,  only  one  matrix-vector  product  is  needed  to  generate  the 
reconstruction.  For  images  of  more  than  modest  sizes,  however,  the  computation  rapidly  becomes  burdensome.  If,  for 
example,  /  and  g  are  100x100,  A  is  a  lO^-by-10  matrix,  and  the  matrix-vector  product  requires  10^  multiplications. 
There  will  also  be  considerable  storage  demands.  Thus,  if  major  computational  resources  are  not  available,  some 
means  of  simplifying  the  calculation  will  be  needed  in  many  cases  of  practical  interest. 

If  the  matrix  A  were  square  circulant  of  order  n(  A  =  ,  subscript  mod  n),  we  could  dramatically  reduce 

the  computational  burden  by  exploiting  the  fact  that  the  Fourier  transform  diagonalizes  a  circulant  If  F  denotes  the 
Fourier  matrix: 


1 

1 

1 

1 

1 

w 

2 

w 

n-1 

w 

Vw 

1 

2 

w 

4 

w 

^2(n-l) 

1 

n-\ 

w 

2(n-l) 
w  ... 

(n-1)^ 

w 

w  =  expiiln/n) 


where 


It  follows  from  Eq.(16)  that 


where 


A  =  F^AF 

A  =  diag  . 


Ap  =  diag  .... 


and  Xj- is  the  complex  conjugate  of  A. „  We  note  that  the  singular  values  [o.]  and  the  eigenvalues  {X^.}  are  related  by 

"rW-  (30) 


The  operational  counts  for  the  SVD  and  the  FFT  are  Oin)  and  O  (nlogn) ,  respectively. 
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5.  STRUCTURE  OF  THE  IMAGING  MATRIX 


Under  some  circumstances  the  imaging  matrix  can  be  readily  modified  to  circulant  form.  For  a  shift-invariant  one¬ 
dimensional  system,  the  image  is  a  convolution  of  the  point  spread  function  and  the  object.  If  the  sampling  intervals  in 
image  and  reconstruction  spaces  are  equal,  the  matrix  A  is  then  Toeplitz  (A=  [a.  _  ^.] ).  If,  in  addition,  the  image  and 
reconstruction  vectors  are  of  equal  length,  A  can  be  padded  to  circulant  form.  In  two  dimensions,  /  and  g  are  matrices 
and  must  be  mapped  into  vectors.  The  manner  in  which  this  is  performed  will  determine  the  structure  of  A.  For 
column-wise  mapping,  for  instance,  again  rvith  equal  sampling  in  image  and  reconstruction  spaces,  A  becomes  block- 
Toeplitz  with  Toeplitz  blocks.  If  the  image  and  reconstruction  matrices  have  the  same  number  of  elements,  A  can  be 
padded  to  become  block-circulant  with  circulant  blocks,  and  the  problem  is  again  amenable  to  Fourier  transform 
methods.  In  both  one  and  two  dimensions,  it  should  be  noted  that  the  penalty  associated  with  the  expansion  of  A  to 
circulant  or  block-circulant  form  is  a  relaxation  of  the  support  constraint  on  the  reconstruction,  which  renders  the 
calculation,  and  in  particular  the  degree  of  resolution  enhancement  achieved,  much  more  sensitive  to  noise  An 
alternative  construction  in  the  two-dimensional  case  is  to  zero-pad /and  g  into  larger  matrices  and  then  to  use  a  linear 
congruential  scan  to  map  the  padded  matrices  into  vectors.  A  then  becomes  a  circulant  matrix  since  the  linear 
congruential  scan  is  an  isomorphism  between  2D  convolution  and  ID  convolution 

For  less  structured  imaging  matrices  (e.g.,  if  the  system  is  weakly  space-variant)  it  may  be  asked  whether 
accelerated  computation  of  the  matrix-vector  product  is  still  possible.  In  this  context,  recent  work  on  circulant 
approximations  to  matrices  of  quite  general  form  appears  highly  relevant,  and  includes  the  following  result.  For  any 
matrix,  A^  say,  we  can  write 

a 

^  L(x^)C^(yJ  (31) 

m  =  1 

where  is  a  circulant  matrix  with  the  same  last  row  as  At,Hx^)  is  a  lower  triangular  Toeplitz  matrix  with 
as  its  first  column,  and  C(y  )  is  a  circulant  matrix  whose  last  row  is  y  .The  {x  }  and  {y  ]  may  be  obtained 
from  the  truncated  SVD  of  the  cyclic  displacement  of  Ap : 

a 

m  =  1 

where  E  is  the  cyclic  downshift  matrix 

'o  0  0 ...  0  r 
1  0  0  ...  0  0 

E  =  0  1  0  ...  0  0  •  (33) 

0  0  0  ...  1  0 


For  imaging  matrices  with  strongly  Toeplitz  features,  a  should  be  a  small  number. 
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6.  THE  REGULARIZED  PSEUDOINVERSE  DECONVOLUTION  ALGORITHM 


Consider  a  two-dimensional  optical  imaging  system  whose  point  spread  function  is  both  time-  and  space- 
invariant.  In  the  presence  of  additive  noise,  the  imaging  equation  connecting  the  input  (extended  object),  the  output 
(degraded  image),  and  the  point  spread  function  (impulse  response)  is  given  by  the  following  two-dimensional 
convolutional  integral  equation 


i(x,y) 


oo  oo 


0  {x', y')p{x-x',y-y')  dx'dy'  +  n  (x,  y) . 


(34) 


In  shorthand  notation 

i  =  p**o  +  M  (35) 

where  o(x,y)  represents  the  extended  object,  i(x,y)  the  degraded  image,  p{x,y)  the  point  spread  function,  and  n(x,y)  the 
noise.  It  is  assumed  that  the  noise  is  independent  of  position  in  the  image. 

The  discrete  version  of  Eq.(34)  can,  of  course,  be  cast  ‘  into  the  following  vector-matrix  form 


g  =  Af+r.  (36) 

When,  in  particular,  /,  g,  and  r  represent  the  one-dimensional  column  vectors  formed  by  stacking  the  rows  or  columns 
of  the  discretized  versioirs  of  the  input  o(x,y),  output  i(x,y),  and  the  noise  n{x,y),  respectively,  the  two-dimensional 
imaging  matrix  A  is  block-Toeplitz  with  Toeplitz  blocks.  It  can  be  shown,  using  known  results  ^  that  the  regularized 
pseudoinverse  Ap  is  block-persymmetric  with  persymmetric  blocks.  The  inverse  of  a  Toeplitz  matrix  is  persymmetric, 
and  persymmetric  matrices  obtained  by  inverting  Toeplitz  matrices  have  much  more  Toeplitz-like  structure  than 
general  persymmetric  matrices.  In  particular,  their  displacement  rank  is  the  same  as  that  of  the  parent-Toeplitz 
matrix^'  Displacement  rank,  it  should  be  noted,  is  a  quantitative  measure  of  the  closeness  of  a  given  matrix  to  being 
Toeplitz.  In  one  dimension,  the  displacement  rank  of  a  Toeplitz  matrix  is  <  2,  and  the  displacement  rank  of  Ap  is  <  4 . 

In  the  early  stages  of  this  investigation,  one  of  the  authors  (R.RB.)  noted  that  for  a  variety  of  point  spread  functions 
being  studied  the  corresponding  computer-generated  regularized  pseudoinverse  matrices  appeared  to  have  banded 
block-Toeplitz  structure.  He  considered  the  implication  of  this  observation.  In  particular,  if  a  space-invariant  point 
spread  function  used  in  a  two-dimensional  linear  convolutional  imaging  equation  leads  to  a  block-Toeplitz  imaging 
matrix  with  Toeplitz  blocks,  then  the  converse  must  also  be  true.  That  is,  given  a  block-Toeplitz  imaging  matrix 
containing  Toeplitz  blocks,  then  the  corresponding  space-invariant  point  spread  function  which  gave  rise  to  this 
imaging  matrix  could  be  easily  ascertained.  This  implies  that  from  the  regularized  pseudoinverse  an  inverse  point 
^read  function  dp  (x,  y)  could  be  constructed  which  could  be  used  to  process  the  image  i(x,y)  and  form  an  estimate 
6p  (x,  y)  of  the  onginal  object.  This  two-dimensional  linear  convolution  technique  is  summarized  by  the  equation 

6p  =  d^-k-ki.  (37) 
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The  technique  was  tested  on  a  digital  computer  with  encouraging  results.  Further  experimental  investigations  indicate 
that  results  obtained  with  this  Regularized  Pseudoinverse  ^convolution  (RAPID)  algorithm  are  comparable  in 
quality  to  those  obtained  using  parametric  Wiener  filtering.  An  error  analysis  is  given  in  Appendix  A.  Results  using 
degraded  images  processed  with  both  the  RAPID  algorithm  and  parametric  Wiener  filtering  are  presented  in  section  8. 

7.  DECONVOLUTION  VIA  WIENER  FILTERING 


Wiener  filtering  is  a  well-known  technique  for  processing  images  degraded  and  corrupted  by  noise  as  described 
by  Eq.(34).  From  a  knowledge  of  the  point  spread  function  p(x,y)  characterizing  the  optical  imaging  system,  it  is 
possible  to  compute  the  corresponding  optical  transfer  function  p(/^./  )  using  the  two-dimensional  Fourier 
transform 


jj 


—00—00 


p(x,y)  exp  [-i2jt  (xf^  +  yf^)  ]  dxdy. 


(38) 


From  the  optical  transfer  function  and  knowledge  of  the  power  spectra  of  object  (/^,/  )  and  noise  if^,f  ) ,  the 
parametric  Wiener  filter  can  be  constructed,  namely 


P(fx’fy) 


|p(4./jP+y 

S„(4./y)  - 

1  'X  'y  1 

- 

J 

(39) 


When  Y  =  1 ,  Eq.(39)  reduces  simply  to  the  Wiener  filter.  If  y  is  variable  we  refer  to  this  as  the  parametric  Wiener  filter. 
In  the  absence  of  noise,  either  form  of  the  Wiener  filter  reduces  to  the  ideal  inverse  filter.  When  the  power  spectra  are 
not  known,  which  is  often  the  case  in  practice,  Eq.(39)  can  be  approximated  by 


P 

<4-V  ^ 

P(4./y)[|P(4./y 

)  ^  +  y] 

(40) 


From  the  Wiener  filter,  using  either  Eqs.(39)  or  (40),  an  inverse  point  spread  function  w^(,x,  y)  can  be  constructed 
using  the  two-dimensional  inverse  Fourier  transform.  That  is. 


v}y(x,  y) 


oo  oo 


^4’ f  +  '2’t(^4  +  y/y)  ]  dj/fy. 


(41) 


With  the  inverse  point  spread  function  given  by  Eq.(41),  an  estimate  dy{x,  y)  of  the  object  can  be  computed  using 
the  two-dimensional  linear  convolutional  equation 

Oy  =  xVy-kifi  (42) 


where,  again,  i(x,y)  is  the  degraded  image. 
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8.  RESULTS 


Preliminary  results  obtained  using  the  RAPID  algorithm  are  presented  in  this  section.  For  purposes  of  comparison, 
results  obtained  using  the  parametric  Wiener  filter  algorithm  are  also  included.  The  optical  system  considered  for  these 
studies  was  a  simple  spherical  converging  lens  as  shown  in  Fig.  1. 


Point  Source 


Figure  1.  A  Simple  Spherical  Converging  Lens  Imaging  System 

The  object  and  image  planes  are  coplanar  and  orthogonal  to  the  optical  axis  of  the  lens.  An  arbitrary  point  source 
located  in  the  object  plane  gives  rise  to  an  intensity  distribution  (the  point  spread  function)  in  the  image  plane.  A 
charge-coupled  device  (CCD),  for  example,  can  be  used  to  measure  the  point  spread  function.  A  ray-trace  program  was 
used  in  the  synthesis  of  four  different  point  spread  functions  dominated  by  spherical  aberration  (Fig.  2),  coma  (Fig.  3), 
astigmatism  (Fig.  4),  and  defocus  (Fig.  5).  The  object  plane  distance  (rrun),  image  plane  distance  (mm),  focal  length 
(mm),  F-number,  tangential  field-angle  (deg.),  and  sagittal  field-angle  (deg.)  associated  with  each  of  these  four  figures 
are  sunrunarized  in  Table  1.  The  same  CCD  model  was  used  in  all  simulations.  The  array  consisted  of  a  31-by-31  planar 
arrangement  of  square  detectors  measuring  0.01  mm  on  a  side  with  a  center-to-center  spacing  of  0.01  mm.  Each  pwint 
spread  function  was  obtained  by  tracing  20,000  rays. 
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Table  1 


Spherical 

Aberration 

Coma 

Astigmatism 

Defocus 

Object  plane  distance 

48.6 

48.6 

48.6 

48.6 

Image  plane  distance 

52.0 

52.0 

50.4 

51.0 

Focal  length 

24.3 

24.3 

24.3 

24.3 

F-number 

4.80 

5.70 

5.70 

4.00 

Tangential  field-angle 

0.00 

2.87 

5.91 

0.00 

Sagittal  field-angle 

0.00 

2.87 

0.00 

0.00 

On  the  top  line,  center  diagram,  of  Figs.  2  through  5  are  mesh  plots  of  the  four  point  spread  functions  considered. 
Each  point  spread  function  p(x,y)  is  represented  by  a  31-by-31  matrix.  The  ideal  extended  object  o(x,y)  used  in  this 
analysis  was  a  256-by-256  matrix  which  is  displayed  as  an  8-bit  gray-level  diagram  in  the  upper-left  hand  comer  of 
each  of  these  figures.  Performing  a  two-dimensional  convolution,  see  Eq.  (35),  of  the  extended  object  with  the  point 
spread  function  (in  the  absence  of  noise)  yields  a  286-by-286  degraded  image  i(x,y)  of  which  the  central  256-by-256 
portion  of  the  degraded  image  is  shown  in  the  upper-right  hand  comer  of  each  of  these  figures.  The  full  286-by-286 
matrix  is  used,  however,  in  subsequent  image  processing  computations. 

On  line  two  of  each  of  the  four  figures  are  mesh  plots  of  the  inverse  point  spread  functions  (31-by-31  matrices) 
d^(x,y)  (left-diagram)  and  w^{x,y)  (right-diagram)  obtained  using  the  RAPID  and  Wiener  filter  algoritluns, 
respectively.  Performing  a  two-dimensional  convolution  of  the  point  spread  function  p(x,y)  with  each  of  the  inverse 
point  spread  functions  d^{x,y)  and  w^(x,y)  yields  processed  point  spread  functions  (61-by-61  matrices).  The 
central  31-by-31  portions  of  these  processed  point  spread  functions  are  shown  as  mesh  plots  on  line  three  (left-  and 
right-diagrams,  respectively).  For  purposes  of  comparison,  a  31-by-31  null-array  with  a  single  non-zero  entry  for  the 
center  pixel  (Delta)  is  also  displayed  on  line  three,  center  diagram. 


The  values  of  the  regularization  parameters  (P  and  y)  which  gave  rise  to  the  best  processed  pioint  s|>read  functions, 
judged  visually,  for  this  analysis  are;  ^herical  aterration  (OandO),  coma  (3x10  '^andlxlO  ^),  astigmatism 
(2xl0~^and2xl0~^),  and  defocus  (5x10  and  1x10  ).  These  same  regularization  parameter  values  were  used  in 

computing  the  object  estimates  dp(x,y)  and  dyix,y)  using  Eqs.  (37)  and  (42),  respectively.  These  object  estimates 
(316-by-316  processed  images)  are  displayed  as  ^ay-level  diagrams  on  the  bottom-line  of  each  of  the  four  figures. 
Only  the  central  256-by-256  portions  of  the  processed  images  are  shown. 


The  results  presented  in  Figs.  2  through  5  were  based  on  studies  using  synthesized  point  spread  functions.  We 
were  fortunate  to  obtain  real  digitized  degraded  images  of  the  planet  Saturn  taken  with  the  wide-field  planetary 
camera  of  the  Hubble  Space  Telescope.  The  upper  diagram  in  Fig.  6  shows  a  400-by-250  degraded  (unprocessed)  image 
of  the  planet  Saturn.  The  second  line  in  Fig.  6  shows  two  31-by-31  degraded  images  of  different  stars  (called  star  #1  and 
star  #2)  also  taken  with  the  wide-field  planetary  camera.  The  RAPID  algorithm  was  used  to  process  the  degraded 
image  of  Saturn  using  the  two  star  images  as  point  spread  functions  characterizing  the  degradation  process.  In 
particular,  star  #1  image  was  used  as  the  input  point  spread  function.  Both  star  #1  and  star  #2  images  were  first 
processed  using  the  inverse  point  spread  function  obtained  using  the  star  #1  image  only.  The  regularization  parameter, 
P,  selected  was  the  one  which  gave  rise  to  processed  star  images  of  equal  quality,  based  on  a  minimum  entropy 
criterion.  This  same  inverse  point  spread  function  was  then  used,  via  Eq.(35),  to  process  the  degraded  Saturn  image. 
The  size  of  the  recoirstructed  image  was  430-by-280.  The  central  4{X)-by-250  portion  of  this  reconstruction  is  shown  in 
the  lower  diagram  of  Fig.  6. 
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Figure  2.  Point  Spread  Function  Dominated  by  Spherical  Aberration 
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Processed  PSF  -  Pseudoinverse  Delta  Processed  PSF  -  Wiener 


Figure  3.  Point  Spread  Function  Dominated  by  Coma 
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Processed  PSF  -  Pseudoinverse 


Delta 


Processed  PSF  -  Wiener 


Processed  -  Pseudoinverse 


Processed  -  Wiener 


Figure  4.  Point  Spread  Function  Dominated  by  Astigmatism 
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Processed  PSF  -  Pseudoinverse 


Delta 


Processed  PSF  -  Wiener 


Processed  -  Pseudoinverse 


Processed  -  Wiener 


Figure  5.  Point  Spread  Function  Dominated  by  Defocus 
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Unprocessed  Saturn  Image 


Star  #1  Image  Star  #2  Image 

■  ■  .  ' 

Processed  Saturn  Image 

Figure  6.  Hubble  Space  Telescope  Imagery 
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9.  CONCLUSIONS 


A  new  algorithm,  the  Regularized  Pseudoinverse  Deconvolution  (RAPID)  algorithm,  has  been  developed  and  has 
xjen  shown  to  be  equivalent  in  performance  to  Wiener  filtering  for  shift-invariant  point  spread  functions.  The 
ilgorithm  can  be  implemented  either  by  direct  convolution  or  indirectly  using  Fast  Fourier  Transform  (FFT) 
echniques.  The  advantage  of  this  approach  is  that,  through  the  use  of  linear  congruential  scanning  and  the  application 
jf  the  circulant  expansion  of  equation  (31),  regularized  reconstruction  methods  can  be  applied  to  extended  images 
degraded  by  shift-variant  point  spread  functions. 
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APPENDIX  A  -  TRADEOFFS  BETWEEN  ACCURACY  AND  SPEED 


In  this  appendix  we  examine  the  RAPID  reconstruction  scheme  in  the  continuous  domain.  An  expression  for  the 
error  displays  a  compromise  between  accuracy  and  speed.  Let  A  denote  the  convolution  operator  on  l"  (31^) : 


Af{x)  =  (a-kj)  (X)  =  J  a(x-y)f(y)dm(}/). 

Let  denote  the  two-dimensional  Fourier  transform  and  let  M.  denote  the  multiplication  operator:  M.cp  =  d(p.  As 
is  well-known  (see  reference  A1  below),  the  Fourier  transform  diagonalizes  the  convolution  operator: 
^  ^  ~  ,with  II  All  =  11  all  ^  .In  this  formalism,  the  regularization  operator  satisfies: 

A -1 „+  A  ,, 

—  M  2  _i  -  ■ 

The  RAPID  operator  is  obtained  by  determining  the  regularization  pseudoinverse  and  then  restricting  it  to  a  region 
near  the  origin.  Introduce  the  associated  convolution  operator 


Ayfix)  =  J*  (w^a)  (x-y)f(y)dmiy), 

9(2 

where  denotes  the  boxcar  window:  (jc)  =  1  if  x  €  [-b,  b]  x  [-b,  b] ,  zero  otherwise.  The  regularization 
operator  for  Ap  is 

A-1  .-H  A  _  w 


Then  the  error  between  the  regularization  operator  and  the  RAPID  operator  is 


a  ^  ^ 

— - - w.-k — - - 

|a|^  +  P  ^  +  P  ^ 

As  b  tends  to  infinity,  the  error  tends  to  zero.  However,  a  large  value  for  b  corresponds  to  an  increased  computational 
load.  Thus,  a  trade  off  must  be  made  between  speed  and  accuracy.  In  addition,  the  error  bound  also  indicates  that  a 
different  window  (such  as  a  square  Hamming)  might  improve  the  accuracy  of  the  regularization  scheme. 


j  +  -  + 

^P"\p 


Al.  N.  Young,  An  Introduction  to  Hilbert  Space,  Cambridge  University  Press,  1988. 
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ABSTRACT 

In  a  paper  presented  at  the  first  coirference  in  this  series,  the  problem  was  considered  of  reconstructing  an  object 
from  image  data  degraded  by  a  compact  linear  operator.  Simulated  synthetic  aperture  radar  measurements  were  used 
to  demonstrate  that,  using  Tikhonov  regularization  techniques,  resolution  well  beyond  the  conventional  limit  can  be 
achieved.  These  image  restoration  methods  have  been  further  extended  and  are  being  investigated  for  application  to 
space-based  surveillance  systems  at  optical  wavelengths  and  also  to  laser  radar  imaging  and  detection  problems.  The 
computational  burden  associated  with  the  reconstruction  procedure  is  of  critical  importance,  since  the  calculation  must 
be  performed  rapidly  enough  to  permit  an  effective  response.  For  realistically-sized  images,  however,  it  can  easily  be 
shown  that  the  operations  coimt  is  prohibitively  lai^,  and  some  method  of  accelerating  the  deconvolution  process 
must  be  found. 

Approximations  to  die  reconstruction  matrix  based  on  drculant  expansions,  which  permit  fast  Fourier  transform 
techniques  to  be  exploited,  represent  a  promising  area  for  research.  A  new  algorithm  for  fast  restoration  of  extended 
images,  the  Regularized  Pseudoinverse  Deconvolution  (RAPID)  algorithm,  suggested  by  the  Toeplitz-like  structure  of 
the  reconstruction  matrix,  is  also  proposed  and  discussed.  Applications  to  both  optical  and  radar  s)^tems  are 
considered  and  illustrated.  In  the  radar  case,  the  squared  modulus  of  the  transmitted  signal's  ambiguity  function  plays 
the  role  of  the  optical  point  spread  function. 

1.  INTRODUCTION 

In  the  absence  of  any  other  degrading  effects,  the  performance  of  an  optical  system  is  ultimately  restricted  by 
diffraction.  The  finite  extent  of  the  entrance  pupil  imposes  a  fundamental  upper  limit  on  the  system's  spatial  frequency 
response.  The  image  quality  of  most  operational  systems  will  not,  however,  approach  this  theoretic.-’  limit  very  closely. 
It  is  possible  that  the  design  or  construction  will  be  flawed,  as  in  the  case  of  the  Hubble  Space  Telescope,  through 
dcfcciivc  manufacture,  assembly  or  quality  assurance  procedures.  The  methods  of  image  restoration  considered  here 
were  originally  aimed  at  extending  the  performance  of  an  optical  system  beyond  the  diffraction  limit^  but  are  in  fiict 
capable  of  compensating  also  for  aberrations  induced  by  the  optical  components.  They  are  inherently  robust  and 
amenable  to  parallel  implementations. 


The  assumption  is  made  that  the  imaging  system  can  be  characterized  by  a  point  spread  function  which  induces  a 
possibly  time-  and  shift- variant  blurring  of  the  image.  It  will  be  assumed  that  the  point  spread  function  is  known  or 
can  be  measured.  In  the  case  of  the  Hubble  wide-field  planetary  camera,  for  example,  there  arc  two  primary 
components  to  the  point  spread  function;  one  due  to  diffraction  by  the  aperture  obstructions  and  the  other  due  to 
spherical  aberration.  For  a  delay  Doppler  imaging  radar,  on  the  other  hand,  the  point  spread  function  is  a  bilinear 
functional  of  the  transmitted  signal. 

The  overall  imaging  system  equation  can  be  written  as  a  Fredholm  integral  equation  of  the  first  kind.  Solutions  to 
ill-posed  problems  of  this  type  are  known  to  be  numerically  unstable^'  Additionally,  it  is  anticipated  that  the  image 
will  be  spatially  sampled  which  will  introduce  discretization  and  associated  noise  processes.  Thus,  after  scanning  the 
image  into  a  vector,  the  integral  representing  the  continuous  image  can  be  rewritten  as  a  matrix  expression.  In  general, 
the  presence  of  the  sensor  noise  takes  the  measured  image  vector  out  of  the  span  of  the  columns  of  the  kernel  matrix, 
which  is  typically  highly  ill-conditioned.  Thus,  when  it  is  desired  to  estimate  what  the  image  would  have  been  in  the 
absence  of  aberrations  or  with  less  diffraction  than  the  instrument  can  provide,  techniques  derived  from  regularization 
theory  are  required  to  restore  stability  to  the  reconstmrtion.  By  introducing  a  suitable  error  criterion  (based  on,  e.g.,  a 
vector  norm),  images  can  be  constructed  which  are,  in  terms  of  the  chosen  criterion,  closer  to  the  undistorted 
geometrical  image  of  the  object  than  the  detected  image  data. 

To  the  extent  permitted  by  the  noise  jr.  the  image,  in-band  effects  can  usually  be  removed  by  some  form  of 
pseudoinverse  filter‘d.  However,  detector  pixellation  and  the  finite  aperture  of  any  system  set  fundamental  resolution 
limits  in  the  performance  of  such  lilters,  and  a  method  for  achieving  spectral  extrapolation  has  to  be  devised.  The 
spatial  spectrum  of  the  object  is  the  Fourier  transform  of  its  amplitude,  in  the  coherent  case,  or  its  intensity,  in  the 
incoherent  case.  If  the  object  is  known  to  be  of  finite  extent,  its  Fourier  transform  is  an  analytic  function,  and  the  out-of- 
band  part  of  the  spectrum  can  in  principle  be  fully  recovered  by  analytic  continuation^  of  the  image  spectrum  after 
removal  of  any  in-band  distortion.  The  inverse  Fourier  transform  of  this  extended  spectrum  would  then  yield  a  perfect 
image  of  the  original  object.  Equivalently,  one  could  attempt  to  solve  directly  the  equation  describing  the  imaging 
process.  This,  however,  involves  the  inversion  of  an  ill-conditioned  matrix,  and  the  restoration  process  is  intrinsically 
unstable,  even  small  amounts  of  noise  rendering  the  results  meaningless.  These  difficulties  may  be  surmounted  by 
applying  the  methods  of  regularization  theory^,  developed  to  deal  with  ill-posed  problems  of  this  type;  the  solution  is 
derived  by  means  of  a  constrained  least-squares  procedure  in  which  a  regularization  parameter  plays  an  essential  role. 
Stability  in  the  restored  image,  which  is  computed  via  the  regularized  pseudoinverse  of  the  imaging  operator,  is 
controlled  by  this  parameter,  whose  optimal  value  depends  on  the  signal-to-noise  ratio  in  the  data. 

2.  NATURE  OF  THE  PROBLEM 

We  wish  to  estimate  the  object  /  from  an  image  g,  given  that 

g  =  Af+r  .(1) 

where  A  is  the  imaging  operator  and  r  represents  the  corrupting  effect  of  additive  noise.  For  clarity  in  the  analysis, 
we  consider  the  one-dimensional  case  with  the  operator  A  given  in  integral  form  by 

(A/)  (y)  =  J  A(x,y)f(x)dx,  c<y<i/.  (2) 

In  the  absence  of  noise,  Eq.(l)  becomes  a  Fredholm  equation  of  the  first  kind,  in  which  the  unknown  function  appears 
only  under  the  integral  sign.  We  can  identify  the  sources  of  difficulty  in  solving  this  equation  in  the  presence  of  noise  or 
other  perturbations  (such  as  computer  rnund-off  error)  by  means  of  a  singular  function  analysis ' . 


We  expand  the  kernel  of  the  integral  in  terms  of  the  singular  functions  (x)  and  v.  (y) ,  orthonnrmal  systems  in 


object  and  image  spaces  respectively,  and  the  singular  values  a.: 


Aix,y)  =  u.(x)  t;.(y). 


The  object  and  image  can  be  expanded  in  the  forms 


«,U) 

i  -  1 

OO 

giy)  =  '^gi  v.{y) 


where  the  coefficients  are  related  to /(x)  and  g(y)  by  the  integral  formulae 

/,•  =  I  f(x)u.(x)dx 

and  J  ^ 

^,=  1  g(y)v.(y)dy. 

In  the  noiseless  case  ( r  =  0) ,  we  find  from  Eqs.  (2)  and  (6) 


{Af)  (y)  =  /.  i;.(y). 


It  follows,  using  Eqs.(5)  and  (8),  that 


and  hence 


g-  =  f- 

I ’t 


5-.  V 


Thus  the  object  function  can  in  principle  be  perfectly  reconstructed  from  the  set  }  of  image  coefficients. 

Now  consider  the  effects  of  noise.  By  expanding  r(y)  in  terms  of  the  v-  (y),  we  can  derive  the  contribution  of  the 
noise  to  the  new  image  coefficients: 


Theeslimate  /(x)  is  now 


g’ .  =  CT.  /.  +  r.. 
*  j  t  I 


J(X)  =  fix)  +  ^  ^ 


(12) 


Image  formation  is  characteristically  described  by  an  integral  transform  of  convolution  type,  i.e.,  >4  is  a 
convolution  operator.  Its  singular-value  spectrum  typically  decays  asymptotically  at  an  exponential  rate  Since  the  r. 
will  in  general  decrease  less  quickly,  the  sum  in  Eq.(12)  will  be  divergent  and  no  bound  will  exist  for  the  'distance'  (in 
the  sense  of  some  appropriate  metric)  between  the  true  object  and  the  reconstruction.  The  effect  of  the  noise  on  the 
reconstructed  image  is  a  manifestation  of  the  fact  that  convolution  is  a  strongly  smoothing  process  -  closely  similar 
images  can  correspond  to  widely  differing  objects.  Thus  image  restoration  is  an  ill-posed  problem,  small  perturbations 
in  the  data  causing  large  changes  in  the  solution  represented  by  Eq.(12). 


3.  THE  REGULARIZED  SOLUTION 


Methods  of  regularization  theory^'  ^  can  be  exploited  to  convert  the  problem  to  a  related  well-posed  one,  i.e.,  one 
for  which  the  solution  exists,  is  unique  and  depends  continuously  on  the  data.  Since  we  shall  later  be  concerned  with 
computed  reconstructions,  we  henceforth  consider  the  problem  in  its  finite  discretized  form;  thus  the  imaging  operator 
A  becomes  a  matrix.  Although  strictly  we  should  introduce  new  symbols,  for  convenience  we  continue  to  use/,  g,  and 
A  to  denote  the  discrete  forms  of  object,  image,  and  imaging  operator.  Generally  and  A  e  C”’ 

To  regularize  the  problem,  we  shall  modify  A.  We  impose  constraints^*^  on  possible  solutions  /'  by  requiring  that 

||A/'-g||2<e2  (13) 

where  e  is  some  suitable  measure  of  the  noise  in  the  image,  and  that 

II  (14) 

where  q  is  some  suitable  measure  of  the  permitted  'signal  strength'  of  the  solution.  ( |1  •  ||  denotes  norm  in  the 
Hilbert  spaces  associated  with  object  and  image.)  We  combine  these  constraints  and  minimize 


ll^/’-5ll^  +  PII  /'  11^ 

where  the  regularization  parameter  p  is  given  by 

P  =  .  (15) 

The  minimum-norm  solution  to  this  constrained  least-squares  problem  is  given  by 


/p  - 

where 

Ap  =  (A^A  +  p/)'^A^. 


(16) 

(17) 


We  note  the  relationship  of  Ap 
pscudoinverse'l 


(which  we  shall  call  the  regularized  pseudoinverso)  to 

A^  =  Urn  An. 

P->0  1^ 


A* ,  the  Moore- Penrose 
(IS) 


H  H 

'The  inverse  of  (A^A  +  p/)  always  exists,  since  A  A  is  non-negative  definite  and  the  regularization  parameter  P 
is  positive.  A  value  of  P  should  be  chosen  which  balances  data  fidelity  against  smoothness  in  the  a'construction. 
Methods  are  also  available  for  determining  P  from  the  image  data  themselves’^' 


4.  SINGULAR  VALUE  AND  FOURIER  DECOMPOSITIONS 


It  is  often  convenient  to  compute  the  regularized  pseudoinverse  via  the  singular  value  decomposition  (SVD)  of  A: 


where^** 


A  -  UZr 


U^U  =  I  ,  V^V  =  VV^  =  I 
tn  n 


E  =  diag(cs^,a^,...,a^),  a.>0. 


The  singular  values  are  assumed  to  have  been  arranged  in  descending  order  of  magnitude 


a,  >CT.  >a-  >  ...  >CT  . 
12  3  n 


Then  we  find 


where 


/p  =  1/  1+  U»  X 


Z^-dag...^  . 

I  «,•  +  P 

This  representation  is  useful  when  the  behavior  of  the  reconstruction  as  a  function  of  the  regularization  parameter  is 
being  studied.  Regularization  can  equivalently  be  achieved  by  setting  p  to  zero,  and  simply  truncating  the  singular 
value  series  at  a  point  which  is  dependent  on  the  noise  leveP^. 

An  ab  initio  computation  of  /p  via  Eq.(23)  requires  the  SVD  of  A  followed  by  two  matrix-vector  multiplications.  If 
the  regularized  pseudoinverse  can  be  precomputed,  only  one  matrix-vector  product  is  needed  to  generate  the 
reconstruction.  For  images  of  more  than  modest  sizes,  however,  the  computation  rapidly  becomes  burdensome.  If,  for 
example,  /  and  g  are  lOO-by-lOO,  A  is  a  lO^-by-10^  matrix,  and  the  matrix-vector  product  requires  10^ 
multiplications.  There  will  also  be  considerable  storage  demands.  Thus,  if  major  computational  resources  are  not 
available,  some  means  of  simplifying  the  calculation  will  be  needed  in  many  cases  of  practical  interest. 

If  the  matrix  A  were  square  circulant  of  order  n(  A  =  ,  subscript  mod  n),  we  could  dramatically  reduce 

the  computational  burden  by  exploiting  the  fact  that  the  Fourier  transform  diagonalizes  a  circulant'^.  If  F  denotes  the 
Fourier  matrix: 

r  1  1  I  ...  1  1 


pH  _  I 

F  -  -j=  1  w 

Vn 


2(«-l) 


w  =  cxp(i2n/n) 


,  n-1  2(n-l) 

1  w  zv 


Under  some  circumstances  the  imaging  matrix  can  be  readily  modified  to  circulant  form.  For  a  shift-invariant  one¬ 
dimensional  system,  the  image  is  a  convolution  of  the  point  spread  function  and  the  object.  If  the  sampling  intervals  in 
image  and  reconstruction  spaces  are  equal,  the  matrix  A  is  then  Toeplitz  (y4  =  [Uj  _  )  •  If/  in  addition,  the  image  and 
reconstruction  vectors  are  of  equal  length,  A  can  be  padded  to  circulant  form.  In  two  dimensions,/  and  g  are  matrices 
and  must  be  mapped  into  vectors.  The  manner  in  which  this  is  performed  will  determine  the  structure  of  A.  For 
column-wise  mapping,  for  instance,  again  with  equal  sampling  in  image  and  reconstruction  spaces,  A  becomes  block- 
Toeplitz  with  Toeplitz  blocks.  If  the  image  and  reconstruction  matrices  have  the  same  number  of  elements,  A  can  be 
padded  to  become  block-circulant  with  circulant  blocks,  and  the  problem  is  again  amenable  to  Fourier  transform 
methods.  In  both  one  and  two  dimensions,  it  should  be  noted  that  the  penalty  associated  with  the  expansion  of  A  to 
circulant  or  block-circulant  form  is  a  relaxation  of  the  support  constraint  on  the  reconstruction,  which  renders  the 
calculation,  and  in  particular  the  degree  of  resolution  enhancement  achieved,  much  more  sensitive  to  noise  .  An 
alternative  construction  in  the  two-dimensional  case  is  to  zero-pad  /and  g  into  larger  matrices  and  then  to  use  a  linear 
congruential  scan  to  map  the  padded  matrices  into  vectors.  A  then  becomes  a  circulant  matrix  since  the  linear 
congruential  scan  is  an  isomorphism  between  2D  convolution  and  ID  convolution^^. 

For  less  structured  imaging  matrices  (eg.,  if  the  system  is  weakly  space-variant)  it  may  be  asked  whether 

accelerated  computation  of  the  matrix-vector  product  is  still  possible.  In  this  context,  recent  work  on  circulant 

approximations  to  matrices  of  quite  general  form'^  appears  highly  relevant,  and  includes  the  following  result.  For  any 

matrix,  /lit  say,  we  can  write 

P  o 

~  ^  ^  ) 

p  0  m 

m  =  1 


(31) 


where  C  is  a  circulant  matrix  with  the  same  last  row  as  An,Hx  )  is  a  lower  triangular  Toeplitz  matrix  with  x 
as  its  first  column,  and  C  is  a  circulant  matrix  whose  last  row  is  y^^^.Thc  and  {y^^j}  may  be  obtained 

from  the  truncated  SVD  of  the  cyclic  displacement  of  /\p  : 


a 

'Lvl 

m  =  1 


where  E  is  the  cyclic  downshift  matrix 


£ 


0  0  0  ...  0  1 

1  0  0  ...  0  0 

0  1  0  ...  0  0  • 


0  0  0  ...  1  0 


(32) 


(33) 


For  imaging  matrices  with  strongly  Toeplitz  features,  a  should  be  a  small  number. 

6.  THE  REGULARIZED  PSEUDOINVERSE  DECONVOLUTION  ALGORITHM 


Consider  a  two-dimensional  optical  imaging  system  whose  point  spread  function  is  both  time-  and  space- 
invariant.  In  the  presence  of  additive  noise,  the  imaging  equation  connecting  the  input  (extended  object),  the  output 
(degraded  image),  and  the  point  spread  function  (impulse  response)  is  given  by  the  following  two-dimensional 
convolutional  integral  equation 


if.x,  y) 


oo  oo 


0  {x',  y')p(x-x\y-  y’)  dx' dy'  +  n{x,y). 


(34) 


In  shorthand  notation 


i  =  pif-ko  +  n  (35) 

where  o(x,y)  represents  the  extended  object,  i(x,y)  the  degraded  image,  p(x,y)  the  point  spread  function,  and  n(x,y)  the 
noise.  It  is  assumed  that  the  noise  is  independent  of  position  in  the  image. 

The  discrete  version  of  Eq.(34)  can,  of  course,  be  cast^®  into  the  following  vector-matrix  form 


g  =  Af+r. 


(36) 


When,  in  particular,  /,  g,  and  r  represent  the  one-dimensional  column  vectors  formed  by  stacking  the  roivs  or  columns 
of  the  discretized  versions  of  the  input  oix,y),  output  i(x,y),  and  the  noise  n(x,y),  respectively,  the  two-dimensional 
imaging  matrix  A  is  block-Toeplitz  with  Toeplitz  blocks.  It  can  be  shown,  using  known  results  that  the  regularized 
pseudoinverse  Ap  is  block-persymmetric  with  persymmetric  blocks.  The  inverse  of  a  Toeplitz  matrix  is  persymmctric, 
and  persymmetric  matrices  obtained  by  inverting  Toeplitz  matrices  have  much  more  Toeplitz-like  structure  than 
general  persymmetric  matrices.  In  particular,  their  displacement  rank  is  the  same  as  that  of  the  parent-Toeplitz 
matrix‘s'  .  Displacement  rank,  it  should  bo  noted,  is  a  quantitative  measure  of  the  closeness  of  a  given  matrix  to  being 
Toeplitz.  In  one  dimension,  the  displacement  rank  of  a  Toeplitz  matrix  is  <  2,  and  the  displacement  rank  of  /Ip  is  <4. 

In  the  early  stages  of  this  investigation,  it  was  noted  that  for  a  variety  of  point  spread  functions  being  studied  the 
corresponding  computer-generated  regularized  pseudoinverse  matrices  appeared  to  have  banded  block-ToepIitz 
structure.  The  implications  of  this  observation  were  considered.  In  particular,  if  a  space-invariant  point  spread  function 
used  in  a  two-dimensional  linear  convolutional  imaging  equation  leads  to  a  block-Toeplitz  imaging  matrix  with 
Toeplitz  blocks,  then  the  converse  must  also  be  true.  That  is,  given  a  block-Toeplitz  imaging  matrix  containing  Toeplitz 
blocks,  then  the  corresponding  space-invariant  point  spread  function  which  gave  rise  to  this  imaging  matrix  could  be 
easily  ascertained.  This  implies  that  from  the  regularized  pseudoinverse  an  inverse  point  spread  function  rfp(x,  y) 
could  be  constructed  which  could  be  used  to  process  the  image  i(x,y)  and  form  an  estimate  op  (x,  y)  of  the  original 
object.  This  two-dimensional  linear  convolution  technique  is  summarized  by  the  equation 

«p  =  ^p*^f.  (37) 

The  technique  was  tested  on  a  digital  computer  with  encouraging  results.  Further  experimental  investigations  indicate 
that  results  obtained  with  this  Regularized  Pseudoinverse  Deconvolution  (RAPID)  algorithm  are  comparable  in 
quality  to  those  obtained  using  parametric  Wiener  filtering.  Results  using  degraded  images  processed  vvitli  both  the 
RAPID  algorithm  and  parametric  Wiener  filtering  are  presented  in  section  8. 

7.  DECONVOLUTION  VIA  WIENER  FILTERING 


Wiener  filtering^^  is  a  well-known  technique  for  processing  images  degraded  and  corrupted  by  noise  as  described 
by  Eq.(34).  From  a  knowledge  of  the  point  spread  function  p(x,y)  characterizing  the  imaging  system,  it  is  possible  to 
compute  the  corresponding  transfer  function  p(f  ,f  )  using  the  two-dimensional  Fourier  transform 

^  y 


P(fx’fy)  = 


—  00—00 


p(x,y)  exp  [-iln  {xf  +  yj  )  ]  dxAy. 

^  y 


(38) 


From  the  transfer  function  and  knowledge  of  the  power  spectra  of  object  ) 

parametric  Wiener  filter  can  be  constructed,  namely 


and  noise  S  (f  ) ,  the 
n  Tv  hf 


When  y  =  1 ,  Eq.{39)  reduces  simply  to  the  Wiener  filter.  If  y  is  variable  we  refer  to  this  ns  the  parametric  Wiener  filter. 
In  the  absence  of  noise,  either  form  of  the  Wiener  filter  reduces  to  the  ideal  inverse  filter.  When  the  power  spectra  are 
not  known,  which  is  often  the  case  in  practice,  Eq.(39)  can  be  approximated  by 


^y(4>/J  = 


p(4./)[1p(4.4)1'+y] 


(40) 


From  the  Wiener  filter,  using  either  Eqs.(39)  or  (40),  an  inverse  point  spread  function  w^(x,y)  can  be  constructed 
using  the  two-dimensional  inverse  Fourier  transform.  That  is. 


II 


— oo — oo 


^y  ^4’  V  ^  ^  ^4  ^  ^  ‘^4'^4 


(41) 


With  the  inverse  point  spread  function  given  by  Eq.(41),  an  estimate  dy(x,  y)  of  the  object  can  be  computed  using  the 
two-dimensional  linear  convolutional  equation 


4  = 


(42) 


u'here,  again,  i(x,y)  is  the  degraded  image. 


8.  APPLICATIONS  TO  THE  OPTICAL  CASE 


Preliminary  results  obtained  for  optical  systems  using  the  RAPID  algorithm  have  been  previously  presented^**  and 
are  summarized  here.  For  purposes  of  comparison,  some  results  obtained  using  the  parametric  Wiener  filter  algorithm 
are  also  included.  The  optical  system  considered  for  these  studies  was  a  simple  spherical  converging  lens.  Tlie  object 
and  image  planes  are  parallel  and  orthogonal  to  the  optical  axis  of  the  lens.  An  arbitrary  point  source  located  in  the 
object  plane  gives  rise  to  an  intensity  distribution  (the  point  spread  function)  in  the  image  plane.  A  charge-coupled 
device  (CCD),  for  example,  can  be  used  to  measure  the  point  spread  function.  A  ray-trace  program  was  used  in  the 
synthesis  of  two  different  point  spread  functions  dominated  by  astigmatism  and  defocus.  The  object  plane  distance 
(mm),  image  plane  distance  (mm),  focal  length  (mm),  F-number,  tangential  field-angle  (deg.),  and  sagittal  field-angle 
(deg.)  associated  with  these  point  spread  functions  are  summarized  in  Tabic  1.  The  same  CCD  model  was  used  in  both 
simulations.  The  array  consisted  of  a  3Tby-31  planar  arrangement  of  square  detectors  measuring  0.01  mm  on  a  side 
with  a  ccnter-to-center  spacing  of  0.01  mm.  Each  point  spread  function  was  obtained  by  tracing  20,(XX)  rays. 

Table  1 


Astigmatism 

Defocus 

Object  plane  distance 

48.6 

48.6 

Image  plane  distance 

50.4 

51.0 

Focal  length 

24.3 

24.3 

F-number 

5.70 

4.00 

Tangential  field-angle 

5.91 

().(K) 

Sagittal  field-angle 

0.00 

().(M) 

On  the  top  line,  center  diagram,  of  Figs.l  and  2  are  mesh  plots  of  the  two  of  the  point  spread  functions  considered. 
Each  point  spread  function  p(x,y)  is  represented  by  a  31-by-31  matrix.  The  ideal  extended  object  o(x,y)  used  in  this 
analysis  was  a  256-by-256  matrix  which  is  displayed  as  an  8-bit  gray-level  diagram  in  the  upper-left  hand  corner  of 
each  of  these  figures.  Performing  a  two-dimensional  convolution,  see  Eq.(35),  of  the  extended  object  with  the  point 
spread  function  (in  the  absence  of  noise)  yields  a  286-by-286  degraded  image  i(x,y)  of  which  the  central  256-by-256 
portion  of  the  degraded  image  is  shown  in  the  upper-right  hand  corner  of  each  of  these  figures.  The  full  286-by-286 
matrix  is  used,  however,  in  subsequent  image  processing  computations. 

On  line  two  of  each  of  those  figures  are  mesh  plots  of  the  inverse  point  spread  functions  (31-by-31  matrices) 
ifpfa:, y)  (left-diagram)  and  w^{x,y)  (right-diagram)  obtained  using  the  RAPID  and  Wiener  filter  algorithms, 
respectively.  Performing  a  two-dimensional  convolution  of  the  point  spread  function  p(x,y)  with  each  of  the  inv'erse 
point  spread  functions  dp(x, y)  and  w^(x,y)  yields  processed  point  spread  functions  (61-by-61  matrices).  The 
central  31-by-31  portions  of  these  processed  point  spread  functions  arc  shown  as  mesh  plots  on  line  throe  (left-  and 
right-diagrams,  respectively).  For  purposes  of  comparison,  a  31-by-31  null-array  with  a  single  non-zero  entry  for  the 
center  pixel  (Delta)  is  also  displayed  on  line  three,  center  diagram. 

The  values  of  the  regularization  parameters,  (3  and  7,  which  gave  rise  to  the  best  processed  point  spread  functions, 
judged  visually,  for  this  analysis  are:  astigmatism,  2x10  ^and2xl0~^,  and  defocus,  SxlO^^and  1x10  These  same 
regularization  parameter  values  were  used  in  computing  the  object  estimates  ()p(x,  y)  and  dy(x,y)  using  Eqs.  (37) 
and  (42),  respectively.  These  object  estimates  (316-by-316  processed  images)  arc  displayed  as  gray-level  diagrams  on 
the  bottom-line  of  each  of  these  figures.  Only  the  central  256-by-256  portions  of  the  processed  images  arc  shown. 

The  results  presented  in  Figs.  1  and  2  were  based  on  studies  using  synthesized  point  spread  functions.  We  were 
fortunate  to  obtain  real  digitized  degraded  images  of  the  planet  Saturn  taken  with  the  wide-field  planetary  camera  of 
the  Hubble  Space  Telescope.  The  upper  diagram  in  Fig.  3  shows  a  4()0-by-25()  degraded  (unprocessed)  image  of  the 
planet  Saturn.  Here  the  point  spread  function  is  weakly  space-variant.  Approximate  allowance  can  bo  made  for  its 
variation  across  the  image  of  Saturn  by  using  the  images  of  two  stars  in  the  same  part  of  the  field -of- view,  obtained 
with  the  telescope  on  another  occasion.  The  second  line  in  Fig.  3  shows  the  31-by-31  degraded  images  of  these  two 
stars  (called  star  #1  and  star  #2).  The  RAPID  algorithm  was  used  to  process  the  degraded  image  of  Saturn  using  the  star 
images  as  point  spread  functions  characterizing  the  degradation  process.  In  particular,  star  #1  image  was  used  as  the 
input  point  spread  function.  Both  star  #1  and  star  #2  images  were  first  processed  using  the  inverse  point  spread 
function  obtained  using  the  star  #1  image  only.  The  regularization  parameter,  p,  selected  was  the  one  which  gave  rise 
to  processed  star  images  of  equal  quality,  based  on  a  minimum  entropy  criterion.  The  entropy  of  the  normalized  image 
is  -Up  where  p^.  is  the  content  of  the  i,j  th  pixel.  This  same  inverse  point  spread  function  was  then  used,  via 

Eq.(37),  to  process  the  degraded  Saturn  image.  The  size  of  the  reconstructed  image  was  43()-by-28().  The  central  4()0-by- 
250  portion  of  this  reconstruction  is  shown  in  the  lower  diagram  of  Fig.3. 

9.  RADAR  IMAGE  RECONSTRUCTION 

Microwave  radar  signal  processing  is  a  well-established  field'S  in  which  the  role  of  the  matched  filter  is  pre- 
i-mincnt.  More  recently,  microwave  imaging  radars  have  become  important  for  terrain  mapping  not  only  of  the  earllr*’ 
but  also  of  other  planetary  objecls^^.  Most  terrain-mapping  radars  are  of  the  synthetic  aperture  type;  the  ground  is 
assumed  to  be  stationary  and  the  radar  is  assumed  to  move  relative  to  it  with  a  known  motion,  so  that  the  radar  return 
can  be  processed  simultaneously  in  range  and  Doppler.  After  pulse  compression,  ratlar  range  (sl.iut  range)  cati  he 
com. erted  to  range  along  the  ground  using  the  assumed  orientation  of  the  r.ular  with  respect  to  the  grouiHl.  After 


Doppler  processing,  cross-range  can  be  computed  from  the  processed  signals.  In  a  simple  side-looking  synthetic 
aperture  radar,  such  as  is  used  in  commercial  aircraft  ground  mapping,  range  and  Doppler  processing  may  be 
performed  independently  and  the  point  spread  function  of  the  imaging  radar  treated  as  separable.  Under  these 
conditions,  image  formation  can  be  viewed  as  a  Kronecker  product  of  two  one-dimensional  processes.  In  a  paper 
presented  at  the  first  conference  in  this  serics^^'  an  image  restoration  technique  for  a  synthetic  aperture  radar  system, 
based  on  the  regularized  pseudoinverse  of  the  imaging  operator,  was  considered. 

In  this  paper,  we  discuss  an  inverse  synthetic  aperture  laser  radar  example,  for  which  the  target  will  be  modeled  as 
a  doubly-spread  (in  range  and  Doppler)  extended  object.  We  will  initially  assume  that  a  monostatic  laser  radar  system 
is  looking  into  space  at  a  distribution  of  reflecting  spheres  of  various  sizes  at  different  ranges  and  associated  with 
different  Doppler  frequencies.  It  is  reasonable  to  assume  also  that  the  target  of  the  laser  radar  is  rougli  relative  to  the 
illuminating  radiation  wavelength.  This  assumption  is  approximately  true  even  for  microwave  radars  operating  in  a 
planetary-mapping  mode“^.  Following  Van  Trees'  standard  statistical  radar/sonar  model'^,  the  return  is  assumed  to  be 
a  Gaussian  random  variable.  The  detection  statistic  used  is  the  log-likelihood  function  and  is  equal  to  the  magnitude 
si]uared  output  of  a  bandpass-matched  filter  or  correlator. 

For  those  unfamiliar  with  the  mathematical  formalism  describing  the  properties  and  performance  of  a  radar 
system,  we  introduce  some  terminology.  The  traditional  approach  to  modeling  the  radar  image  is  via  the  asymmetric 
ambiguity  function'^: 

OO 

X(v.'t)  =  fz*  (f)z  (f  +  T)c.xp(  +  i27tv/)i//  (43) 


where  2(f)  =  s(f)  +i.s(f)  and  s(f)  denotes  the  Flilbert  transform  of  s(/)  and  x  and  v  represent  delay  and 
Doppler  shift.  This  narrowband  ambiguity  function  may  be  viewed  as  the  response  of  either  a  matched  filter  receiver 
to  a  point  scatterer  as  a  function  of  the  target  delay  and  Doppler  shift,  or  a  bank  of  matched  filters  to  a  point  scatterer  at 
one  particular  delay  and  Doppler  shift.  The  symmetric  form  of  this  ambiguity  function  is 


OO 


/4(v,x) 


cxp(  +  iKVX)x(,v,z) 


X  X 


+  -  )cxp(  +  i2K\'t)dt. 


—  OO 


(44) 


Feig  and  Griinbaum^^'  have  made  use  of  the  symmetric  ambiguity  function  to  demonstrate  the  connection  between 
radar  detection  and  tomography.  Bernfeld^^  first  recognized  the  analogy  between  the  ambiguity  function  (more 
precisely,  its  squared  modulus)  and  certain  mathematical  entities  appearing  in  the  equations  governing  x-ray 
tomographic  data.  Snyder  et  al?^  noted  the  analogy  with  lime-of-llight  positron  emission  tomography.  |See  also  ref. 
33.|  Whitehouse  and  Boashaslr^'*  established  a  related  analogy  betwevn  the  radar  image  and  that  obtained  by  time-of- 
flight  positron  emission  tomography  for  the  bistatic  case,  where  the  phase  of  the  transmitted  signal  is  not  needed. 


The  Wigner-Ville  distribution  (WVD)  of  a  real  signal  s(t),  associated  with  the  complex  analytic  signal  z(t),  is  a  time- 
frequency  distribution  defined  as'^'^ 

OO 


W(f,/) 
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zit  +  ^)z*  (t  exp  ( -;2  nfx  )flx. 


(4.^) 


—  OO 


Mote  that  the  WVD  is  the  double  Fourier  transform  of  the  symmetric  ambiguity  function.  The  WVD  is  always  real, 
u’hereas  the  ambiguity  function  is  in  general  complex.  For  a  given  signal,  the  WVD  is  related  to  the  squared  modulus 
of  the  ambiguity  function  by  a  double  convolution^^: 


W**W  =  |>li  .  (46) 

As  stated  above,  the  radar  return  is  assumed  to  be  a  Gaussian  random  variable,  and  the  detection  statistic  used  is  the 
log-likcliho('d  function.  Under  this  stochastic  treatment  of  the  radar  rctlection,  the  determination  of  the  target's 
scattering  function  a  becomes  a  problem  in  deconvolution: 


E(|X„j2)  =  |Z„ 


a 


where  £  is  the  expectation  operator,  s  denotes  the  transmitted  signal,  r  denotes  the  received  signal  and  x  's  the  cross¬ 
ambiguity  function  between  transmitted  and  received  signals: 


X^^(v,T)  =  j  s*  (t)r(t  +  x)exp(  +  i2itvt)iit.  (48) 

— OO 

Thus  we  can  see  from  Eq.(47)  that  the  known  function  plays  the  role  ('f  the  point  spread  function  in  the  radar 

imaging  case.  The  radar  deconvolution  problem  can  also  be  stated  in  the  closely-related  torm 

£(Wp  =  W^-k-ka  (49) 

where  and  arc  the  WVDs  of  the  transmitted  and  received  signals  respectively.  A  consistent  estimator  for  this 
equation  is  discussed  in  Appendix  A. 

We  observe  that  the  total  energy,  %  in  the  ambiguity  function  is 


Jj 


'E  =  I  I  |/4(v,  x)\  dvdx 


and  that  |  A]  possesses  the  self-transform  property^^: 


\A(\/,x)\^exp[i2n(xf-vt)]dvdx  =  (A(/,/)|“. 


We  also  note  an  important  distinction  between  the  optical  and  radar  cases.  The  resolution  of  a  corrected  optical 
system  of  modest  F-number  (ratio  of  focal  length  to  objective  diameter)  can  be  improved  simply  by  increasing  the 
aperture.  Suppose  an  iris  is  illuminated  by  a  point  source  whose  radiant  intensity  can  be  adjusted  to  maintain  the 
power  passing  through  the  aperture  at  a  constant  level.  As  the  diameter  of  the  iris  increases,  the  central  peak  of  the 
point  spread  function  becomes  both  higher  and  narrower;  the  height  increases  in  proportion  to  the  area  of  the 
aperture  while  the  width  decreases  in  proportion  to  its  diameter.  Note  that  the  contribution  of  the  volume  of  the  central 
peak  relative  to  the  integrated  sidelobes  remains  constant.  For  the  radar  svsiem,  on  the  other  hand,  it  can  be  shown 
that,  for  a  constant-energy  signal  of  increasing  time-bandwidth  product,  the  width  of  the  central  peak  in  the  ambiguity 
function  decreases,  while  its  magnitude  remains  constant  and  the  integrated  volume  under  the  sidelobes  increases. 


Tlius,  radar  systems  become  increasingly  aberrated  as  the  time-bandwidth  product  increases,  and  we  are  obliged  to 
resort  to  post-processing  of  the  data  to  maintain  feature  contrast  as  vve  increase  the  resolution  for  distributed  objects. 
This  situation  brings  to  mind  the  Hubble  wide-field  planetary  camera,  whose  design  resolution  was  achieved,  but 
whose  performance  was  degraded  by  the  enlargement  of  the  point  spread  function  sidelobes  through  severe  spherical 
aberration. 

10.  APPLICATION  TO  THE  RADAR  CASE 

As  an  illustration  of  the  application  of  the  ideas  discussed  here  to  the  radar  case,  we  consider  the  ambiguity 
function  associated  with  a  laser  radar  system^®  radiating  pulses  with  Gaussian  envelopes.  The  ambiguity  function  is 
then  a  two-dimensional  Gaussian.  Since  the  WVD  of  the  signal  is  also  a  two-dimensional  Gaussian,  the  computation 
involved  in  the  image  restoration  procedure  has  exactly  the  same  form  for  both  F:qs.(47)  and  (49).  The  computer 
software  used  in  generating  the  simulation  results  in  Figs.  1  and  2  was  applied  to  the  case  of  a  circularly  symmetric 
Gaussian  point  spread  function.  The  results  of  this  simulation  are  summarized  in  Fig.4. 
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Figure  2.  Point  Spread  Function  Dominated  by  Defocus 
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Figure  3.  Hubble  Space  Telescope  Imagery 
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Figure  4.  Gaussian  Point  Spread  Function 


APPENDIX  A  -  A  CONSISTENT  ESTIMATION  FOR  VVVD  RADAR  IMAGING 

Consider  the  radar  deconvolution  problem  represented  by  Eq.(49).  This  equation  assumes  the  expected  value  (4' 
;he  WVD  c'f  the  recei\ed  signal  is  known.  This  appendix  discusses  consistent  estimation  Cor  this  equation  using  a 
pulse-averaging  scheme,  but  depends  on  substantial  statistical  assumptions  regarding  the  random  field  of  radar 
reflectors. 

Under  the  assumptions  outlined  in  Section  9,  the  received  signal  {  r  (f)  }  is  a  stochastically  filtered  version  ot  the 
transmitted  signal  {s(t)}  (Eq.  (8)  of  [36]): 


r(t) 


oo  oo 


D  (x,y)s  (t  -  X)  exp  i+iZnyt]  dxdij. 


—  CO  —  oo 


(Al) 


Here  D(x,y)  denotes  the  random  field  determined  by  the  reflectivity  coefficient  of  the  scatterers  as  a  function  of  radial 
range  x  and  radial  Doppler  y.  For  ease  of  exposition,  vve  are  assuming  both  {rtf)}  and  { s  t  f)  }  are  analytic  and  have 
subsumed  a  deterministic  phase  term  in  D(x,y).  Then  the  WVD  for  the  received  signal  has  the  form: 


W(t,f)  = 


OO  CO  OO  OO  OO 


y^ )  D*  (x^,  y^)s(t  +  --  .V, )  .s»  ^  -  .v^) 
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The  reflectivity  field  is  assumed  to  be  independent  in  range  and  Doppler: 


£[D(x^,  D*  (X2.y2)J  =  a (x^,  y^) 5(x^  -  X2)  6  fy^  -y2) . 


(A2) 


(A3) 


By  taking  the  expectation  of  Eci-  (A2)  and  using  the  covariance  of  D(x,y)  in  Eq.  (A3),  it  is  straight-forward  to  obtain 
Eq.(49); 


E  I  l\M/,  f)  I 


OO  00 


a  ( X,  y )  W^  ( f  -  .V.  f  -  y)  dxd  1/ 
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It  is  one  thing  to  write  down  an  expectation  and  quite  another  to  have  a  consistent  estimator.  The  remainder  of  this 
appendix  considers  a  pulse-averaging  estimator  for  Eq.  (49).  Suppose  N  time-vlelayed  pulses  are  transmitted  so  that  the 
nth  received  waveform  is  given  by 


r 

n 


(/) 


oO  OO 


D^(x,y)cxpl  +  i2nyt}s^ 


(t  -  x)dxdy, 
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where  the  nf/i  pulse  is  s^(t)  =  s  (f  -  f^)  and  the  associated  random  fields  D^(x,y)  satisfy: 


I 

I 

(D-1)  Each  (x,  y)  is  a  zero-mean,  circular  Gaussian  random  field. 

(D-2)  The  sequence  [D^(x,y)}  is  independent  and  identically  distributed  with  respect  to  n. 

(D-3)  Each  D^(x,  y)  exhibits  zero  correlation  with  respect  to  a:  and  y. 

Assumptions  (D-1),  (D-2),  and  (D-3)  imply  that  Eq.(A3)  holds  for  each  (x,  y)  for  a  single  scattering  function 
a  (that  is,  a  =  cr^  )  while  the  expectation  between  fields  is  zero:  £  [D^(X|,yp  (x,,,  y2)]  =0  (or  m^n. 

Consider  averaging  of  time-delayed  WVD  of  these  received  pulses: 

N 

^  Y  (»  +  <„./)•  (A6) 

n  =  \ 

It  is  straight-forward  to  verify  that  the  mean  of  this  estimator  is 

E[W(N-,t,f)]  =  (A7) 

The  covariance  computation  is  more  involved  but  it  can  be  bounded: 


^CovlW(N;lyf^),  W(N;/2,/2)l|5^l)5)l2X(|a||J. 


(A8) 


2 

Thus,  W(N;t,f)  is  a  consistent  estimator  for  Eq.  (49)  as  N  provided  the  transmitted  signal  belongs  to  L  (91)  and  the 
scattering  function  a  belongs  to  L^(9?x9?).  While  W(N;t,p  is  consistent,  note  that  assumption  (D-2)  begins  to  break 
down  for  large  N.  Indeed,  while  each  D  (x,  y)  may  still  be  Gaussian  and  independent  with  respect  to  n,  it  seems 
unlikely  that  the  associated  scattering  functions  could  be  identical  as  It  is  conjectured,  for  a  frequency 

independent  reflectivity  coefficient,  that  an  alternative  consistent  estimator  can  be  obtained  by  averaging  the  WVDs  of 
an  ensemble  of  frequency  shifted  pulses  in  a  manner  similar  to  Eq.(A6).  This  frequency  averaged  estimator  should 
overcome  the  difficulty  of  accounting  for  the  non-stationarity  of  the  stochastic  reflectivity  field  due  to  the  motion  in 
range  of  the  scatterers  during  the  observation  interval  of  the  cn,-<'mblc. 


